{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Understanding the PyBaMM pipeline and performance profile of the Simulation class\n",
    "\n",
    "This notebook discusses the main steps that PyBaMM performs to go from an initial model to a final solution, and how you can use this information to optimise the performance of your simulations. We will also discuss the `pybamm.Simulation` class, which is a convenience class that wraps up the entire pipeline into a single class, and how you can use this class most efficiently to run your simulations, taking into account the performance characteristics of the different steps of the pipeline and the caching behaviour of the `pybamm.Simulation` class.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import numpy as np\n",
    "\n",
    "import pybamm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## The PyBaMM Pipeline\n",
    "\n",
    "The PyBaMM pipeline is a series of steps that are taken to go from a model to a final solution. \n",
    "At the start of the pipeline, a pybamm model consists of a set of equations make up of symbols arranged in an [expression tree](https://docs.pybamm.org/en/stable/source/examples/notebooks/expression_tree/expression-tree.html). Each symbol could be a parameter (`pybamm.Parameter`), a variable (`pybamm.Variable`), a binary operator (`pybamm.Addition`, `pybamm.Multiplication`, etc.), or a unary operator (`pybamm.Negate`, `pybamm.Gradient`). You can see all the possible symbols in the API docs [here](https://docs.pybamm.org/en/stable/source/api/expression_tree/index.html). The pipeline is a series of steps that first transforms this expression tree into a form that can be solved by a numerical solver, and then solves the equations to get the solution.\n",
    "\n",
    "\n",
    "The pipeline is as follows:\n",
    "\n",
    "1. **Parameter replacement** - The first step of the pipeline is performed by the [`pybamm.ParameterValues`](https://docs.pybamm.org/en/stable/source/api/parameters/parameter_values.html#pybamm.ParameterValues) class, which (a) replaces all `pybamm.Parameter` symbols with their corresponding values (e.g. using `pybamm.Scalar`), and (b) replaces all `pybamm.FunctionParameter` symbols with either a `pybamm.Scalar`, `pybamm.Interpolant`, or `pybamm.Function` object. \n",
    "2. **Discretisation** - The next step is to discretise the model equations. This is done by the [`pybamm.Discretisation`](https://docs.pybamm.org/en/stable/source/api/spatial_methods/discretisation.html) class, which replaces all the symbols representing spatial derivatives (e.g. `pybamm.Gradient`) with a matrix that can be used to solve the equations numerically.\n",
    "3. **Solver setup** - The next step is to set up the numerical solver. This is done by the [`pybamm.BaseSolver`](https://docs.pybamm.org/en/stable/source/api/solvers/base_solver.html#pybamm.BaseSolver) class, which transforms the expression trees into a form that can be evaluated by the numerical solver. PyBaMM has a few different forms that can be used, but generally, this means that each expression tree is transformed into a CasADi function that can be evaluated at a given time and state.\n",
    "4. **Solver solve** - The next step is to solve the equations. This is done by the numerical solver, which in most cases is one of the solvers provided by the [Sundials](https://sundials.readthedocs.io/en/latest/#) suite of solvers. The solver takes the CasADi functions generated in the previous step and uses them to solve the equations numerically.\n",
    "5. **Post-processing** - The output of the previous step is the solution to the equations at a set of time points. The final step is to post-process this solution to get the variables of interest. This is done by the [`pybamm.Solution`](https://docs.pybamm.org/en/stable/source/api/solvers/solution.html) and [`pybamm.ParameterValues`](https://docs.pybamm.org/en/stable/source/api/parameters/parameter_values.html#pybamm.ParameterValues) classes, which take the solution and the parameter values and return the variables of interest interpolated at the desired time points.\n",
    "\n",
    "For every simulation you run with PyBaMM, these steps are always performed in order to generate each solution, even if they are hidden from you by various convenience classes like [`pybamm.Simulation`](https://docs.pybamm.org/en/stable/source/api/simulation.html) . Understanding the pipeline is important because it can help you identify where the performance bottlenecks are in your simulations, and how you can optimise them to get the best performance.\n",
    "\n",
    "Below is a simple example that demonstrates the different parts of the pipeline, using the DFN model. While you probably have seen something similar to this script before, the one difference is that we have split the solver setup (step 3) and solve (step 4) parts of the pipeline into two separate steps. Often these are combined into the first call to the solver, but we have split them here to make it clear that they are separate steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLEklEQVR4nO3dd3zT1f4/8FfSkY40adO96aJAoYzKKFspaFUE13VdCw4U3PPrrQPFe7Uo93pFr/byE72gXsSLghurIkVANi2UVWhp6aCldCXpStvk8/ujNBBpQ0rTfpL09Xw88pAmJ8n7kF7yumd8jkQQBAFEREREDkIqdgFERERE1sRwQ0RERA6F4YaIiIgcCsMNERERORSGGyIiInIoDDdERETkUBhuiIiIyKEw3BAREZFDYbghIiIih8JwQ0RERA5F1HCTmZmJxMREKBQKKBQKJCcnY+PGjWaf8/bbbyM+Ph7u7u4IDw/Hk08+iZaWln6qmIiIiGyds5hvHhYWhqVLlyIuLg6CIGD16tWYM2cOcnJykJCQcFH7NWvW4C9/+Qs++ugjTJw4EcePH8f8+fMhkUjw1ltvidADIiIisjUSWzs4U6VSYdmyZbjvvvsueuyRRx7B0aNHsWnTJuN9Tz/9NHbt2oVt27ZZ9PoGgwGnT5+Gl5cXJBKJ1eomIiKiviMIArRaLUJCQiCVmp94EnXk5kJ6vR7r1q1DY2MjkpOTu2wzceJEfPrpp9i9ezfGjRuHkydP4ocffsDdd9/d7evqdDrodDrjz+Xl5Rg2bJjV6yciIqK+V1pairCwMLNtRA83eXl5SE5ORktLC+RyOTZs2NBt+LjzzjtRXV2NyZMnQxAEtLe3Y+HChXj++ee7ff2MjAwsWbLkovtLS0uhUCis1g8iIiLqOxqNBuHh4fDy8rpkW9GnpVpbW1FSUgK1Wo0vvvgCK1euxJYtW7oMONnZ2bj99tvxt7/9DePHj0dBQQEef/xxLFiwAC+99FKXr//HkZvOvxy1Ws1wQ0REZCc0Gg2USqVF39+ih5s/SklJQUxMDFasWHHRY1OmTMGECROwbNky432ffvopHnjgATQ0NFxyDg7o2V8OERER2YaefH/b3HVuDAaDyUjLhZqami4KME5OTgA6FhoRERERibrmJj09HampqYiIiIBWq8WaNWuQnZ2NrKwsAEBaWhpCQ0ORkZEBAJg9ezbeeustjB492jgt9dJLL2H27NnGkENEREQDm6jhpqqqCmlpaaioqIBSqURiYiKysrIwc+ZMAEBJSYnJSM2LL74IiUSCF198EeXl5fD398fs2bPx2muvidUFIiIisjE2t+amr3HNDRERkf2x6zU3RERERL3BcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNxY0en6Zhyr1IhdBhER0YDGcGMlG/MqMG3ZZqSvz+M5V0RERCJiuLGSpEE+kEokyCmpx/aCGrHLISIiGrAYbqwkwMsNd46PAAAs33ScozdEREQiYbixooXTYuDqLMWe4jrsPFkrdjlEREQDEsONFQUq3HDbFeEAgHc2nRC5GiIiooGJ4cbKFk6PgYuTBDtO1mBPMUdviIiI+hvDjZWFervjliSO3hAREYmF4aYPPDQ9Bs5SCbaeqMb+kjqxyyEiIhpQGG76QLjKAzeNCQUAvMvRGyIion7FcNNHHpoeC6kE2Jx/FluOnxW7HCIiogGD4aaPDPLzxC1JYQCA+1fvwfr9ZSJXRERENDAw3PShV+cMx3UjgtGmF/DU/w7g7V94cT8iIqK+xnDTh9xcnPDuHaOxcFoMAODtX07g6XUH0NpuELkyIiIix8Vw08ekUgn+kjoEr984Ak5SCdbvL8cdH+xE4dkGsUsjIiJySAw3/eTO8RH4aP5YyGXO2HeqDqnLt+LdTSc4ikNERGRlDDf9aNpgf2x8fAqmDvZHa7sB//j5OGa/u43XwiEiIrIihpt+Fq7ywOp7xuLt20ZB5emK/DNa3Jz5O5774iCqNC1il0dERGT3GG5EIJFIMHd0KDY9NQ03jwmDIACf7y3F9L9nY/kvJ9DU2i52iURERHZLIgywvckajQZKpRJqtRoKhULscgAA+07V4m/fH0VOST0AIFAhwzOz4nHTmDA4SSXiFkdERGQDevL9zXBjIwRBwHcHK/DGj8dQVtcMAIgNkOOZWYNxdUIQJBKGHCIiGrgYbsyw1XDTqaVNj9W/F+P97EKom9sAAIlhSjx7dTwmx/ox5BAR0YDEcGOGrYebTpqWNnzw20l8uK0ITa16AMC4KBWemjkYE6J9Ra6OiIiofzHcmGEv4aZTdYMO728uxKc7T6FV33FNnORoXzw5czDGRalEro6IiKh/MNyYYW/hplOFuhnvby7E2j0laNN3fGSTY/3w6FWxGM+RHCIicnAMN2bYa7jpVF7fjPc2F2Dd3lJjyBkfpcJjM+IwMcaXa3KIiMghMdyYYe/hplNZXRMyswuxbm+ZcbpqTIQ3Hr0qDtPj/RlyiIjIoTDcmOEo4aZThboZK7acxGe7S6A7d07V0GAFHr4yBqnDg3mdHCIicggMN2Y4WrjpVKVpwQdbT+K/u0qMu6ui/DyxcFo05o4OhczZSeQKiYiILh/DjRmOGm461TW2YtXvxVj1e7HxOjmBChnunRSFO8dHwMvNReQKiYiIeo7hxgxHDzedGnTtWLPrFD7cVoQzGh0AwEvmjLsmROLeSYMQoHATuUIiIiLLMdyYMVDCTSddux5f557Gii2FKDzbCABwcZLghpGhWDA1CkOCHP/vgIiI7B/DjRkDLdx0MhgEbDpWhX9vKcS+U3XG+6fE+eH+KdGYGsejHYiIyHYx3JgxUMPNhfaX1GHl1pP48VAlDOc+/bgAOe6dHIW5o0Lh7srFx0REZFsYbsxguDmvtLYJH24rwrq9pWg8t8PK28MFd46LwN3JkQhWuotcIRERUQeGGzMYbi6maWnD//aUYtXvxSirawYAOEkluCYhCGnJkRgXpeKUFRERiYrhxgyGm+7pDQJ+PnIGH20vwu6iWuP9Q4K8MG/iIMwZFQIPV2cRKyQiooGK4cYMhhvLHK3Q4OMdxdiQU46Wto4rH3u5OePmMWH484RIxAbIRa6QiIgGkp58f0v7qaYuZWZmIjExEQqFAgqFAsnJydi4cWO37adPnw6JRHLR7brrruvHqgeGocEKZNyUiF3pKXjxuqGIUHlA29KOVb8XI+WtLbjj/+3E9wcr0HruyAciIiJbIerIzbfffgsnJyfExcVBEASsXr0ay5YtQ05ODhISEi5qX1tbi9bWVuPPNTU1GDlyJFauXIn58+db9J4cubk8BoOArQXV+HTnKWw6esa4y8pP7opbksJx+9hwDPLzFLdIIiJyWHY9LaVSqbBs2TLcd999l2z79ttvY/HixaioqICnp2VfrAw3vVde34y1u0vw+Z5SVGl1xvsnxfri9rERmJUQyLOsiIjIqnry/W0zq0P1ej3WrVuHxsZGJCcnW/ScDz/8ELfffrvZYKPT6aDTnf8C1mg0va51oAv1dsfTs+Lx2Iw4/HqsCp/tLsGW42exvaAG2wtq4OPhghtHh+G2seGID/ISu1wiIhpgRB+5ycvLQ3JyMlpaWiCXy7FmzRpce+21l3ze7t27MX78eOzatQvjxo3rtt0rr7yCJUuWXHQ/R26sq6yuCZ/vKcW6vWWo1LQY7x8Z7o0/XRGG2SNDoOChnUREdJnsalqqtbUVJSUlUKvV+OKLL7By5Ups2bIFw4YNM/u8Bx98EDt27MDBgwfNtutq5CY8PJzhpo/oDQJ+O34Wn+8pxS9Hz6D93OIcmbMU1wwPwq1J4ZgY4wuplNfNISIiy9lVuPmjlJQUxMTEYMWKFd22aWxsREhICF599VU8/vjjPXp9rrnpP2e1OnyVU451+0px/EyD8f4QpRvmjg7FzUlhiPHnlnIiIro0u1xz08lgMJiMtHRl3bp10Ol0+POf/9xPVdHl8PeSYcHUaNw/JQoHy9T4Yl8Zvs4tx2l1C97PLsT72YUYGe6NW8aE4vrEEPh4uopdMhEROQBRR27S09ORmpqKiIgIaLVarFmzBm+88QaysrIwc+ZMpKWlITQ0FBkZGSbPmzJlCkJDQ7F27doevydHbsTV0qbHpqNVWL+/DNnHz0J/btrKxUmC6fEBuHlMKK4cEsDdVkREZMJuRm6qqqqQlpaGiooKKJVKJCYmGoMNAJSUlEAqNb3OYH5+PrZt24affvpJjJKpl9xcnHBdYjCuSwzGWa0OX+eWY0NOOQ6f1uDnI2fw85EzULq74NoRwZg7KgRjB6m4PoeIiHrE5tbc9DWO3Nim/Eot1ueU4euc0ya7rUKUbpg9KgRzR4ViSJAXD/AkIhqg7HpBcV9juLFteoOAnSdr8HVuOTbmVUKrazc+Fhcgx+yRIbg+MRjRXIhMRDSgMNyYwXBjP1ra9MjOr8JXOafxa36VyTlWw0MVuD4xBNeNCEa4ykPEKomIqD8w3JjBcGOfNC1t+OnwGXx74DS2FVQbFyIDwKhwb1yfGIxrRwQjxNtdxCqJiKivMNyYwXBj/2oadNh4qBLfHTyNXUW1uPA3eHSEN64bEYxrhgchzIcjOkREjoLhxgyGG8dSpW3Bj4cq8d3BCuwpNg06I8OUSB0RjGsSgnhiORGRnWO4MYPhxnGd0bQg63AlfsirwO6iWlwwc4UhQV6YlRCEaxKCMDSYu66IiOwNw40ZDDcDw1mtDlmHK5F1uBI7CmuMZ1wBQLjKHbOGBWHWsEAkRfrA2Ulq5pWIiMgWMNyYwXAz8NQ3tWLT0Sr8eLgSvx0/C90Fu658PFwwY2ggUoYGYkqcHzxlNnciCRERgeHGLIabga2ptR2/Ha/GT0cqseloFdTNbcbHXJ2lmBjjixlDAzFjSAB3XhER2RCGGzMYbqhTu96A3cW1+OVIFX45egYltU0mjw8NVuCqIf64akggRoV7w4nHQBARiYbhxgyGG+qKIAgoqGrAz0fPYNPRKuwvqTPZeaXydMXUOD9cOSQAU+L8oeIJ5kRE/YrhxgyGG7JEbWMrthyvwqajVdhy/Cy0LeePgZBIgJFh3rgyPgDT4v0xIlTJUR0ioj7GcGMGww31VJvegP2n6pB9/Cw2H6vCsUqtyeM+Hi6YEuePaYP9MSXODwEKN5EqJSJyXAw3ZjDcUG9VqJuxJf8sthw/i20nqk0O9wQ6rqkz9VzQGTtIBTcXJ5EqJSJyHAw3ZjDckDW16Q3ILa03hp1Dp9Uma3VkzlKMHaTCpFg/TI71w7AQBaewiIguA8ONGQw31JdqGnTYXliDrcfPYuuJalRqWkwe9/ZwQXK0LybG+mFijC+i/Tx5tWQiIgsw3JjBcEP9RRAEFJ5twLYT1dhWUIOdJ2vQ8IcprGClG5JjfI2BJ5TX1iEi6hLDjRkMNySWdr0BB8rU+L2gGtsLq7H/VD1a9QaTNhEqD0yM8UVyjC/GR/kiSMnFyUREAMONWQw3ZCuaW/XYe6oWOwpr8HthDfLK1dAbTP/nGOXniQnRKkyIZtghooGN4cYMhhuyVdqWNuwp7gg7u4pqcahcjT9kHUT6emB8lArjo3wxPlqFMB8PcYolIupnDDdmMNyQvVA3t2FvcS12nuw+7IR6u2NclArjo1QYF6VCFBcoE5GDYrgxg+GG7JWmpQ37iuuws6gGu07WdjmN5Sd3xdhBHUFn7CAVhgZz6zkROQaGGzMYbshRNOrakVNSj91FNdhZVIvc0nq0tpsuUPaSOSNpkA/GRakwbpAKI8KUkDnzooJEZH8YbsxguCFHpWvX42CZGruLarG7qBb7TtVdtPVc5izFqHBvjI9SYWyUCmMifOApcxapYiIiyzHcmMFwQwNFu96AY5Va7CqqxZ6iWuwurkVtY6tJGyepBMNDlecWKatwxSAVlO4uIlVMRNQ9hhszGG5ooOq8qODuojrsKe4Y3SmvbzZpI5EAQ4MUmBDtiwnRHWt3vD1cRaqYiOg8hhszGG6Iziura8Ke4lrsOtkRdk5WN5o83hl2kmN8MSnWF+OifCHnNBYRiYDhxgyGG6LuVWlasKuoY/v5zpM1KDxrGnacpRKMDPfGpBhfTB3sj1Hh3nB2kopULRENJAw3ZjDcEFmuStuCXSdr8XthDbYXVKOktsnkcYWbM6YM9se0wf6YHu+PAC9eQZmI+gbDjRkMN0SXr7S2CdsLqrG1oBrbTlRD3dxm8vjoCG9cnRCEqxOCEOXnKVKVROSIGG7MYLghso7Og0C35Fch+/hZHCxTmzweFyDHtSOCccOoEMT4y0WqkogcBcONGQw3RH2jUt2Cn4+ewU+HK7GjsAbtF1w9eXioAjeMDMHskSEIVrqLWCUR2SuGGzMYboj6nrq5DZuOnsG3B07jtxPVxmMiJBJgcqwfbkkKw9UJQXBz4dWSicgyDDdmMNwQ9a+aBh1+OFSJb3LLsae4zni/l5szZo8Mwa1JYRgV7s0DP4nILIYbMxhuiMRzqqYRX+4vx5f7ykwuIBgXIMdtY8Mxd3Qo/OQyESskIlvFcGMGww2R+AwGATtP1mDdvjL8kFcB3bkDP52lEswYGoBbksIxPd4fLryGDhGdw3BjBsMNkW3RtLTh2wOn8b89pThwwY4rX09XzBkViluSwjAshP9bJRroGG7MYLghsl3HKjX4cl8ZNuScRnWDznj/kCAvzB4ZgusTgxHpy+vnEA1EDDdmMNwQ2b52vQG/nTiLL/aV4ZcjVWjVG4yPjQxT4vrEEFwzPAjhKg8RqySi/sRwYwbDDZF9UTe1IetwJb49eBrbC6pxweVzMCTICylDAzFzWCBGhCohlXLHFZGjYrgxg+GGyH6d1erw46EKfJ9XgT3Fdcbr5wCAv5cMk2P9MCXOD5Nj/RCg4DlXRI6E4cYMhhsix1Df1IrN+VX4+cgZbMk/i8ZWvcnj8YFemBCtwtgoFcYNUjHsENk5hhszGG6IHI+uXY99p+qw7UQ1tp6oxqHTavzxX7ZBvh5IilRhVIQ3Rod7Iz7Ii1vNiewIw40ZDDdEjq+2sRU7Cmuwu6gGu4vrcKxSc1HYkTlLMTxUiRGhSgwLUSAhRIG4AC+4OjPwENkihhszGG6IBh51cxv2n6rD/pI65JbW40BpPTQt7Re1c3WSIi5QjvggLwwJ8kJ8kAJDgrwQ4CXj8RBEIrObcJOZmYnMzEwUFxcDABISErB48WKkpqZ2+5z6+nq88MILWL9+PWpraxEZGYm3334b1157rUXvyXBDRAaDgKKaRhworcfh0xocPq3G4dMaaLsIPACgdHdBbIAccQFyxJ67xfjLEeLtDifu0CLqF3YTbr799ls4OTkhLi4OgiBg9erVWLZsGXJycpCQkHBR+9bWVkyaNAkBAQF4/vnnERoailOnTsHb2xsjR4606D0ZboioK4IgoLS2GUcqNMiv1CL/jAbHKrUorm402X5+IVdnKaJ8PRHt74koP08M8vPEIF9PDPLzgL+coz1E1mQ34aYrKpUKy5Ytw3333XfRY//+97+xbNkyHDt2DC4uLpf1+gw3RNQTLW16FFU34kRVAwrOaDv+W9WAUzVNJhcX/CNPVyeEqzwQofJApG/Hf8NUHgj38UCYjzvcXJz6sRdE9s8uw41er8e6deswb9485OTkYNiwYRe1ufbaa6FSqeDh4YGvv/4a/v7+uPPOO/Hcc8/Byanrfyh0Oh10uvOXcddoNAgPD2e4IaJe0RsElNc1o7C6ASfPNqK4uhHFNR238rrmbkd7OvnJZQhXuSPE2x1h3u4I9XFHqHfHzyFKdyjcnTnyQ3SBnoQb536qqVt5eXlITk5GS0sL5HI5NmzY0GWwAYCTJ0/i119/xV133YUffvgBBQUFeOihh9DW1oaXX365y+dkZGRgyZIlfdkFIhqAnKQSRPh6IMLXA1fGmz6ma9ejrK4ZJbVNKKlpQkltE07VNKGsrgmltU1obNWjukGH6gYdckrqu3x9D1cnBCvdEOLtjmClG4KUnf91Q7DSDcEKBiCi7og+ctPa2oqSkhKo1Wp88cUXWLlyJbZs2dJlwBk8eDBaWlpQVFRkHKl56623sGzZMlRUVHT5+hy5ISJbIggC1M1tKK1tRmldE07XN6O8vhnldR3/rVC3oLax1aLXcndxQpDSDUEKN5Pg0xmEgpVuUHm6MgCRQ7CrkRtXV1fExsYCAJKSkrBnzx4sX74cK1asuKhtcHAwXFxcTKaghg4disrKSrS2tsLV1fWi58hkMshksr7rABFRD0gkEnh7uMLbwxUjwpRdtmlp06NC3YLT9c04Xd+MSnULKjQtHf9Vt6BC3Yz6pjY0n1sPVFTd2O37uTpLjUEnxPuCqS9vd4R6d9zn4Sr6VwGRVdncb7TBYDAZabnQpEmTsGbNGhgMBkilHRfaOn78OIKDg7sMNkRE9sjNxQlRfh07sLrT0qbHGU1H2KlUt6DyXPg5Xd+MSk0LTte3oLpBh9Z2A07VdEyLdUfl6YrQc8EnzKfjFnpu4XOYjzu83C5vAweRWEQNN+np6UhNTUVERAS0Wi3WrFmD7OxsZGVlAQDS0tIQGhqKjIwMAMCiRYvwr3/9C48//jgeffRRnDhxAq+//joee+wxMbtBRNTv3FycEOnriUjf7gNQa7sBZzQdgadC3YLT6uZzo0EtxukwbUs7ahtbUdvYirxydZevo3R3QbjKHeE+HghXeRh3gUWoPBDq7c6rOpPNETXcVFVVIS0tDRUVFVAqlUhMTERWVhZmzpwJACgpKTGO0ABAeHg4srKy8OSTTyIxMRGhoaF4/PHH8dxzz4nVBSIim+XqLDWGke5oWtpQXteMsrpmlNc1oby+488dtybUNbVB3dwGdXkbDpVrLnq+RAKEKN0RrnJHpMqzY5H1ue3vkb6eULpz1If6n+gLivsbr3NDRGS5Bl37uV1ezSitbULpBX8uqW1Cc5ve7PN9PFzOX9zw3AUOo/3kGOTnweku6hG7vM5Nf2G4ISKyDkEQUN3Q2rHlvbYRp2rOb30vrmlCdUPX6yc7+clliD63tijq3FWeY/zliPT14IntdBGGGzMYboiI+keDrh2najpCT1F1I07VNKK4ugknqxvNBh8nqQSRKg9E+8sRE+CJWP/zZ3pxtGfgYrgxg+GGiEh82pa2c0Gnwbid/eTZRpw824DG1u6nugIVMgwO9MLgQC/EB3ohPsgLcYFybmcfABhuzGC4ISKyXYIg4IxGh8KzDSg823GOV+etStv1aI9EAkSqPDAkSIEhwV4YEqTA0GAvRKg8eAFDB8JwYwbDDRGRfVI3t50LOlocq9Ti+Bkt8isbup3i8nJzxrBgBRJClEgIUSAxTIkYfzmkUgYee8RwYwbDDRGRY6lp0CG/UoujlVocq9DgaKUGxysbujy1XS5zxvBQBUaGeWNUuDeSIn0QoHAToWrqKYYbMxhuiIgcX5vegIKqBhw+rcGhcjUOn1bjULmmy63rYT7uuCLSB0mRPhgX5YvBgXJOZ9kghhszGG6IiAamdr0BBWcbcLBUjdyyeuSU1CO/UgPDH74FfT1dMSHaFxOiVUiO8UOMvyfDjg1guDGD4YaIiDppW9qQW1qPfafqsLe4DvtO1V00uhPq7Y5p8f6YPtgfE2P9IJdxZ5YYGG7MYLghIqLutLYbcLCsHjsKa7DjZA32nqpDa/v5tTsuThJMiPbF1QlBuDohCP5eMhGrHVgYbsxguCEiIks1tbZj18laZOdXIfv4WZPT1SUSYGykCqkjgnB9YgiDTh9juDGD4YaIiC7XybMNyDp8Bj8eqsCBsvOnqDtJJZg+2B83J4VhxtAAyJydRKzSMTHcmMFwQ0RE1lBe34wfD1Xi2wOnkVtab7xf6e6CG0eHYt7EQYjy8xSvQAfDcGMGww0REVlbQVUD1u8vw/r95ajUtADomLa6Kj4A906OwsQYX+646iWGGzMYboiIqK/oDQK2FVTj49+LselYlfH++EAvPHRlDGYnhvAKyZeJ4cYMhhsiIuoPJ882YPXvxVi3rwxN5w4DjQ/0wlOzBmPWsECO5PQQw40ZDDdERNSf1M1t+Pj3Yvy/rSehbWkHACSGKfHMrHhMHewvcnX2g+HGDIYbIiISg7qpDf9vayH+s73YOJJz7YggLL4+AUFKnm91KQw3ZjDcEBGRmKobdHhvcwE+3nEKeoMAucwZT88ajLTkQXDiepxuMdyYwXBDRES24MhpDV74Kg85JfUAgOGhCrxxcyISQpTiFmajevL9Le2nmoiIiOgCw0IU+HLhRLx+4wgo3JxxqFyDG9//HR/vKMYAG3ewOoYbIiIikUilEtw5PgKbnp6OlKGBaG03YPHXh7Ho0/1QN7eJXZ7dYrghIiISmb+XDB+kJWHx9cPg4iTBj4crcd07W5FTUid2aXaJ4YaIiMgGSCQS3Ds5Cl8umogIlQfK6prxpxU78HVuudil2R2GGyIiIhuSGOaN7x6bjGsSgtCmF/D42lys3HpS7LLsCsMNERGRjVG4ueD9u8bg3klRAIC/fX8Ur31/BAYDFxpbguGGiIjIBkmlErx0/VCkpw4BAHywtQhP/S8Xre0GkSuzfQw3RERENkoikeDBaTF4608j4SyV4Kvc03jov/vQrmfAMYfhhoiIyMbdNCYMK+ddAZmzFL8crUL6+jxeC8cMZ0saqVSqHr2oRCLB/v37ERkZeVlFERERkanp8QF4784xePDTfVi3rwz+XjL83zVDxC7LJlkUburr6/H2229Dqbz0JaEFQcBDDz0EvV7f6+KIiIjovJRhgXj9xuF47ss8vJ9diAAvGeafW3RM51kUbgDg9ttvR0BAgEVtH3300csuiIiIiLp329gInNXq8PefjmPJd0fg5yXD9YkhYpdlUywKNwZDzxYuabXayyqGiIiILu3hK2NRpdXh4x2n8OTnuQhWuiMp0kfssmyGxQuKv/vuux6HHCIiIrI+iUSCl2cnGC/098TnOdC28CyqThaHm7lz5yI8PBwvvPACCgoK+rImIiIiugQnqQRv3pqIMB93lNY2Y/HXh8UuyWZYHG6Kiorw4IMPYu3atYiPj8e0adPwySefoLm5uS/rIyIiom4o3Fzw9m2jIJUAG3LKeQ7VORaHm/DwcCxevBiFhYX45ZdfMGjQICxatAjBwcFYuHAh9uzZ05d1EhERUReuGKTCo1fFAQBe3HAIpbVNIlckvsu6iN+VV16J1atXo6KiAsuWLUNeXh4mTJiAkSNHWrs+IiIiuoRHr4pFUqQPtLp2PPF57oC/gnGvrlDs5eWFGTNm4Morr4S3tzeOHDlirbqIiIjIQs5OUrx92yh4yZyx71Qd3ttcKHZJorqscNPc3IyPP/4Y06dPR1xcHNauXYunnnoKxcXFVi6PiIiILBGu8sBf5w4HALy3uWBAT0/1KNzs3LkTDzzwgHGdTVhYGH755RcUFBTghRdeQGhoaF/VSURERJcwZ1QIJsX6olVvwN9/yhe7HNFYHG6GDRuGSZMmYf/+/cjIyEBFRQU+/fRTXHnllX1ZHxEREVlIIpEgPXUoAODr3NM4VK4WuSJxWBxuUlJSsH//fuzduxeLFi2y6JwpIiIi6l/DQ5WYM6rjOIaMjUcH5OnhFoebd955h7uhiIiI7MAzs+Lh6iTF9oIa/HaiWuxy+p1F4WbMmDGoq6uz+EUnT56M8nJeSIiIiEgM4SoPpCVHAgAyfjgKvWFgjd5YdHBmbm4uDhw4AJVKZdGL5ubmQqfTXbJdZmYmMjMzjbusEhISsHjxYqSmpnbZftWqVbjnnntM7pPJZGhpabGoLiIiooHi4Stj8fneUhyr1OKrnHLcnBQmdkn9xqJwAwAzZsyweN5OIpFY1C4sLAxLly5FXFwcBEHA6tWrMWfOHOTk5CAhIaHL5ygUCuTnn18Bbul7ERERDSQ+nq54+MpYLN14DP/4KR/XJQbDzcVJ7LL6hUXhpqioqMcvHBZ26YQ4e/Zsk59fe+01ZGZmYufOnd2GG4lEgqCgoB7XQ0RENNDMnzgIq38vxml1Cz7fU4p5EweJXVK/sCjcREZG9nUd0Ov1WLduHRobG5GcnNxtu4aGBkRGRsJgMGDMmDF4/fXXuw1CAKDT6UymyDQajVXrJiIislVuLk5YND0Gi78+jP/uOoW05MgBMePRq+MXrCEvLw9yuRwymQwLFy7Ehg0bMGzYsC7bxsfH46OPPsLXX3+NTz/9FAaDARMnTkRZWVm3r5+RkQGlUmm8hYeH91VXiIiIbM7c0aFwc5Hi+JkG7C+xfHOQPZMIIm+Ab21tRUlJCdRqNb744gusXLkSW7Zs6TbgXKitrQ1Dhw7FHXfcgb/+9a9dtulq5CY8PBxqtRoKhcJq/SAiIrJVz647gHX7ynDzmDD840/2eVkXjUYDpVJp0fe36CM3rq6uiI2NRVJSEjIyMjBy5EgsX77coue6uLhg9OjRKCgo6LaNTCaDQqEwuREREQ0kd4yPAAB8d/A01E1tIlfT90QPN39kMBgs2kYOdKzTycvLQ3BwcB9XRUREZL9Gh3tjSJAXdO0GbMjpfimHo7iscFNfX4+VK1ciPT0dtbW1AID9+/f3+MJ96enp+O2331BcXIy8vDykp6cjOzsbd911FwAgLS0N6enpxvavvvoqfvrpJ5w8eRL79+/Hn//8Z5w6dQr333//5XSDiIhoQJBIJLjz3OjNZ7tLHf5IBouvc9Pp4MGDSElJgVKpRHFxMRYsWACVSoX169ejpKQEH3/8scWvVVVVhbS0NFRUVECpVCIxMRFZWVmYOXMmAKCkpARS6fn8VVdXhwULFqCyshI+Pj5ISkrC77//btH6HCIiooFszqhQvP7DUeSf0WJ/ST2SIn3ELqnP9HhBcUpKCsaMGYM333wTXl5eOHDgAKKjo/H777/jzjvvNF5t2Fb1ZEESERGRI3lm3QF8sa8MtySF4e+32tfC4j5dULxnzx48+OCDF90fGhqKysrKnr4cERER9ZM7xl2wsLjZcRcW9zjcyGSyLi+Ed/z4cfj7+1ulKCIiIrK+MRHeiA/0QkubAV/lOO4B1z0ONzfccANeffVVtLV1JD6JRIKSkhI899xzuPnmm61eIBEREVmH6cLiEpGr6Ts9Djf/+Mc/0NDQgICAADQ3N2PatGmIjY2Fl5cXXnvttb6okYiIiKxk7uhQuDhJcKxSi6LqRrHL6RM93i2lVCrx888/Y9u2bTh48CAaGhowZswYpKSk9EV9REREZEVKdxdcEanCjpM1yM6vQpRflNglWV2Pw02nyZMnY/LkydashYiIiPrB9Hh/7DhZgy3Hz+KeSQw3eOedd7q8XyKRwM3NDbGxsZg6dSqcnJx6XRwRERFZ3/T4AGRsPIYdhTVoadPDzcWxvrN7HG7++c9/4uzZs2hqaoKPT8cFgOrq6uDh4QG5XI6qqipER0dj8+bNPIGbiIjIBg0OlCNY6YYKdQt2nqzB9PgAsUuyqh4vKH799dcxduxYnDhxAjU1NaipqcHx48cxfvx4LF++HCUlJQgKCsKTTz7ZF/USERFRL0kkEkyP77h8S3b+WZGrsb4eh5sXX3wR//znPxETE2O8LzY2Fn//+9+Rnp6OsLAwvPnmm9i+fbtVCyUiIiLrmTa4I9xsOc5wg4qKCrS3t190f3t7u/EKxSEhIdBqtb2vjoiIiPrEpFg/OEslKKpuxKkax9oS3uNwc+WVV+LBBx9ETk6O8b6cnBwsWrQIV111FQAgLy8PUVGOt/qaiIjIUXi5uRgPz3S0qakeh5sPP/wQKpUKSUlJkMlkkMlkuOKKK6BSqfDhhx8CAORyOf7xj39YvVgiIiKyns6FxNn5VSJXYl09PhW807Fjx3D8+HEAQHx8POLj461aWF/hqeBEREQdjlZokLp8K9xcpMhdPMumt4T35Pv7si/iN2TIEAwZMuRyn05EREQiGxLkhSCFGyo1LdhVVGtcZGzvLivclJWV4ZtvvkFJSQlaW1tNHnvrrbesUhgRERH1LYlEgmmD/fH53lJsyT87cMPNpk2bcMMNNyA6OhrHjh3D8OHDUVxcDEEQMGbMmL6okYiIiPrI9PiOcJN9vAqLMUzscqyixwuK09PT8cwzzyAvLw9ubm748ssvUVpaimnTpuHWW2/tixqJiIioj0yK84OTVIKTZxtRWtskdjlW0eNwc/ToUaSlpQEAnJ2d0dzcDLlcjldffRVvvPGG1QskIiKivqNwc0FSROeWcMfYNdXjcOPp6WlcZxMcHIzCwkLjY9XV1darjIiIiPrF5Dg/AMCe4jqRK7GOHq+5mTBhArZt24ahQ4fi2muvxdNPP428vDysX78eEyZM6IsaiYiIqA+NDPcGAOSVq8UtxEp6HG7eeustNDQ0AACWLFmChoYGfP7554iLi+NOKSIiIjs0IlQJACiqboSmpQ0KNxeRK+qdHoeb6Oho4589PT3x73//26oFERERUf9SeboizMcdZXXNOFSuxsQYP7FL6pUer7mJjo5GTU3NRffX19ebBB8iIiKyH52jN3ll9j811eNwU1xcDL1ef9H9Op0O5eXlVimKiIiI+teIsI5wc9AB1t1YPC31zTffGP+clZUFpVJp/Fmv12PTpk0YNGiQVYsjIiKi/pEY6g0AODSQws3cuXMBdFyqed68eSaPubi4YNCgQTwJnIiIyE51TkudqmmCuqkNSg/7XVRscbgxGAwAgKioKOzZswd+fva92IiIiIjOU3q4INLXA6dqmpBXrjZe+8Ye9XjNTVFREYMNERGRAxoe2rnupl7cQnrJopGbd955x+IXfOyxxy67GCIiIhJPYqgS3x+ssPt1NxaFm3/+858WvZhEImG4ISIislPGHVN2vh3conBTVFTU13UQERGRyDqnpcrqmlHb2AqVp6vIFV2eHq+5uZAgCBAEwVq1EBERkYgUbi6I9vMEYN/nTF1WuPn4448xYsQIuLu7w93dHYmJifjkk0+sXRsRERH1s87RG3ted3NZB2e+9NJLeOSRRzBp0iQAwLZt27Bw4UJUV1fjySeftHqRRERE1D8Sw5T45sBpHCyrF7uUy9bjcPPuu+8iMzMTaWlpxvtuuOEGJCQk4JVXXmG4ISIismOOcMZUj6elKioqMHHixIvunzhxIioqKqxSFBEREYkjIVQJiQQ4rW5BdYNO7HIuS4/DTWxsLP73v/9ddP/nn3+OuLg4qxRFRERE4pDLnBHjLwdgv4uKezwttWTJEtx222347bffjGtutm/fjk2bNnUZeoiIiMi+jAhVoqCqAXllalwZHyB2OT1m8cjNoUOHAAA333wzdu3aBT8/P3z11Vf46quv4Ofnh927d+PGG2/ss0KJiIiof3Suu7HXi/lZPHKTmJiIsWPH4v7778ftt9+OTz/9tC/rIiIiIpEknrtScZ6dnjFl8cjNli1bkJCQgKeffhrBwcGYP38+tm7d2pe1ERERkQiGhSgglQBnNDqc1drfomKLw82UKVPw0UcfoaKiAu+++y6Kioowbdo0DB48GG+88QYqKyv7sk4iIiLqJx6uzgj1cQcAnDzbIHI1Pdfj3VKenp645557sGXLFhw/fhy33nor3nvvPUREROCGG27oixqJiIionw3y7TiGobimUeRKeq5XZ0vFxsbi+eefx4svvggvLy98//331qqLiIiIRBR17oypouomkSvpucsON7/99hvmz5+PoKAgPPvss7jpppuwffv2Hr1GZmYmEhMToVAooFAokJycjI0bN1r03LVr10IikWDu3LmXUT0RERGZ0xluiqvtb+SmR9e5OX36NFatWoVVq1ahoKAAEydOxDvvvIM//elP8PT07PGbh4WFYenSpYiLi4MgCFi9ejXmzJmDnJwcJCQkdPu84uJiPPPMM5gyZUqP35OIiIgubZCf/U5LWRxuUlNT8csvv8DPzw9paWm49957ER8f36s3nz17tsnPr732GjIzM7Fz585uw41er8ddd92FJUuWYOvWraivr+9VDURERHSxqAvW3BgMAqRSicgVWc7icOPi4oIvvvgC119/PZycnKxeiF6vx7p169DY2Ijk5ORu27366qsICAjAfffdZ9FWdJ1OB53u/DY2jUZjlXqJiIgcWZiPO5ylErS0GXBG24JgpbvYJVnM4nDzzTff9EkBeXl5SE5ORktLC+RyOTZs2IBhw4Z12Xbbtm348MMPkZuba/HrZ2RkYMmSJVaqloiIaGBwdpIiXOWBoupGFFU32lW46dVuKWuIj49Hbm4udu3ahUWLFmHevHk4cuTIRe20Wi3uvvtufPDBB/Dz87P49dPT06FWq4230tJSa5ZPRETksAb5egAAiu1sx1SPD860NldXV8TGxgIAkpKSsGfPHixfvhwrVqwwaVdYWIji4mKTdToGgwEA4OzsjPz8fMTExFz0+jKZDDKZrA97QERE5JgG+XkC+WftblGx6OHmjwwGg8kamU5DhgxBXl6eyX0vvvgitFotli9fjvDw8P4qkYiIaEA4f60bhhuLpaenIzU1FREREdBqtVizZg2ys7ORlZUFAEhLS0NoaCgyMjLg5uaG4cOHmzzf29sbAC66n4iIiHrPeJVihhvLVVVVIS0tDRUVFVAqlUhMTERWVhZmzpwJACgpKYFUKvqyICIiogGpc+TmVG2TXW0HlwiCIIhdRH/SaDRQKpVQq9VQKBRil0NERGSz9AYBQ1/6Ea16A7Y9dyXCfDxEq6Un398cFiEiIqIuOUklCFd1bAG3px1TDDdERETULeOiYjvaMcVwQ0RERN2yx0XFDDdERETUrUF2eDo4ww0RERF1i9NSRERE5FA6R25Ka5vQrjeIXI1lGG6IiIioW8EKN8icpWjTCzhd3yJ2ORZhuCEiIqJuSaUSRJ47QNNepqYYboiIiMgse9sxxXBDREREZtnbAZoMN0RERGSWcTs4p6WIiIjIEXBaioiIiBxK57RUaV0z2uxgOzjDDREREZkVqJDB3cUJeoOAsrpmscu5JIYbIiIiMksiuWA7eHWDyNVcGsMNERERXdL5HVNNIldyaQw3REREdEn2dIAmww0RERFdUqSqY1qqrI4jN0REROQAApVuAIAzGp3IlVwaww0RERFdUqBXZ7ix/cMzGW6IiIjokoLOjdzUNLZC164XuRrzGG6IiIjoknw8XODq1BEbzmpte2qK4YaIiIguSSKRIEAhA2D7U1MMN0RERGSRIEXH1FSlmiM3RERE5ADO75jiyA0RERE5AHvZMcVwQ0RERBYJUnasualkuCEiIiJHEKjgyA0RERE5kCCFfVylmOGGiIiILBJo3C3VAkEQRK6meww3REREZJHOqxQ3t+mhaWkXuZruMdwQERGRRdxcnKB0dwEAVNnwuhuGGyIiIrJYoML2d0wx3BAREZHFLlx3Y6sYboiIiMhinTumqmz48EyGGyIiIrIYR26IiIjIoXSeL8U1N0REROQQjNNSDDdERETkCLhbioiIiBxK58jNWa0O7XqDyNV0jeGGiIiILOYrl8FJKoFBAGoaW8Uup0sMN0RERGQxJ6kE/vJzU1M2umOK4YaIiIh6xNZ3TDHcEBERUY8EnVtUbKs7phhuiIiIqEc6FxVz5KYLmZmZSExMhEKhgEKhQHJyMjZu3Nht+/Xr1+OKK66At7c3PD09MWrUKHzyySf9WDEREREFGK9SbJtHMDiL+eZhYWFYunQp4uLiIAgCVq9ejTlz5iAnJwcJCQkXtVepVHjhhRcwZMgQuLq64rvvvsM999yDgIAAXH311SL0gIiIaODpHLk5Y6MjNxJBEASxi7iQSqXCsmXLcN9991nUfsyYMbjuuuvw17/+1aL2Go0GSqUSarUaCoWiN6USERENSNsLqnHXyl2IC5Dj56em9ct79uT722bW3Oj1eqxduxaNjY1ITk6+ZHtBELBp0ybk5+dj6tSp3bbT6XTQaDQmNyIiIrp8tn6VYlGnpQAgLy8PycnJaGlpgVwux4YNGzBs2LBu26vVaoSGhkKn08HJyQnvv/8+Zs6c2W37jIwMLFmypC9KJyIiGpA6TwbXtrSjqbUdHq6ixwkToo/cxMfHIzc3F7t27cKiRYswb948HDlypNv2Xl5eyM3NxZ49e/Daa6/hqaeeQnZ2drft09PToVarjbfS0tI+6AUREdHA4eXmAk9XJwDAGY3tLSq2uTU3KSkpiImJwYoVKyxqf//996O0tBRZWVkWteeaGyIiot676u/ZOFndiM8WTEByjG+fv59drrnpZDAYoNNZngJ72p6IiIh6L9CGd0yJOkmWnp6O1NRUREREQKvVYs2aNcjOzjaOwqSlpSE0NBQZGRkAOtbPXHHFFYiJiYFOp8MPP/yATz75BJmZmWJ2g4iIaMAJUjLcdKmqqgppaWmoqKiAUqlEYmIisrKyjAuES0pKIJWeH1xqbGzEQw89hLKyMri7u2PIkCH49NNPcdttt4nVBSIiogEpwIZ3TNncmpu+xjU3REREvfef7UVY8u0RXDsiCO/fldTn72fXa26IiIjI9p2/SrHtrXtluCEiIqIeO3++lO1NSzHcEBERUY91Liiu0rbAYLCtFS4MN0RERNRjAV4dC4rb9ALqmlpFrsYUww0RERH1mIuTFH5yVwC2t2OK4YaIiIgui5+8Y/SmtpEjN0REROQAfDw6Rm4YboiIiMgh+Hi6AADqm9pErsQUww0RERFdFm+O3BAREZEjUZ0LN/XcLUVERESOwNujY1qqjtNSRERE5AhUnh0jN7zODRERETmEzt1SDDdERETkEIzTUo2cliIiIiIHwGkpIiIiciidW8GbWvVoadOLXM15DDdERER0WRRuznCSSgDY1oX8GG6IiIjoskgkEvgYt4PbztQUww0RERFdts6pqTobukoxww0RERFdNpVxOzinpYiIiMgBdG4Hr+W0FBERETmCzu3g9ZyWIiIiIkfgzWkpIiIiciQqT+6WIiIiIgfibYPnSzHcEBER0WXz4VZwIiIiciTnp6W45oaIiIgcAC/iR0RERA6l8yJ+Wl072vQGkavpwHBDREREl03h7gJJx9mZNnN4JsMNERERXTYnqQRKd9vaDs5wQ0RERL2isrF1Nww3RERE1Cud50tx5IaIiIgcQuf5UrayHZzhhoiIiHrF1q5SzHBDREREveLTOS3FNTdERETkCHw4LUVERESOxNbOl2K4ISIiol7x4ZobIiIiciTGNTecliIiIiJHcH7NDUduiIiIyAF0Tkupm9ugNwgiV8NwQ0RERL3UeYViQegIOGJjuCEiIqJecXGSwsvNGYBtTE0x3BAREVGv2dJ2cFHDTWZmJhITE6FQKKBQKJCcnIyNGzd22/6DDz7AlClT4OPjAx8fH6SkpGD37t39WDERERF1xZYu5CdquAkLC8PSpUuxb98+7N27F1dddRXmzJmDw4cPd9k+Ozsbd9xxBzZv3owdO3YgPDwcs2bNQnl5eT9XTkRERBfysaGTwSWCIIi/rPkCKpUKy5Ytw3333XfJtnq9Hj4+PvjXv/6FtLQ0i15fo9FAqVRCrVZDoVD0tlwiIiIC8OTnudiQU4701CF4cFqM1V+/J9/fzlZ/98uk1+uxbt06NDY2Ijk52aLnNDU1oa2tDSqVqts2Op0OOp3O+LNGo+l1rURERGTq/FWKB/i0FADk5eVBLpdDJpNh4cKF2LBhA4YNG2bRc5977jmEhIQgJSWl2zYZGRlQKpXGW3h4uLVKJyIionNs6WRw0cNNfHw8cnNzsWvXLixatAjz5s3DkSNHLvm8pUuXYu3atdiwYQPc3Ny6bZeeng61Wm28lZaWWrN8IiIigm1dpVj0aSlXV1fExsYCAJKSkrBnzx4sX74cK1as6PY5f//737F06VL88ssvSExMNPv6MpkMMpnMqjUTERGRKVs6PFP0cPNHBoPBZI3MH7355pt47bXXkJWVhSuuuKIfKyMiIqLu2NLhmaKGm/T0dKSmpiIiIgJarRZr1qxBdnY2srKyAABpaWkIDQ1FRkYGAOCNN97A4sWLsWbNGgwaNAiVlZUAALlcDrlcLlo/iIiIBrrOaan6gT5yU1VVhbS0NFRUVECpVCIxMRFZWVmYOXMmAKCkpARS6fllQZmZmWhtbcUtt9xi8jovv/wyXnnllf4snYiIiC5w4W4pQRAgkUhEq0XUcPPhhx+afTw7O9vk5+Li4r4rhoiIiC5b5+GZeoMATUs7lO4uotUi+m4pIiIisn9uLk7wcHUCIP52cIYbIiIisgpb2THFcENERERW4ePZMRVVL/KOKYYbIiIisorOkZtaTksRERGRI+C0FBERETmU8xfyY7ghIiIiB3D+fCmuuSEiIiIHYJyW4pobIiIicgTenJYiIiIiR6LydIWTVAJBELcOmzsVnIiIiOzTpBg/FLyWKuq5UgDDDREREVmJVCpuqOnEaSkiIiJyKAw3RERE5FAYboiIiMihMNwQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNwQERGRQ2G4ISIiIocy4E4FFwQBAKDRaESuhIiIiCzV+b3d+T1uzoALN1qtFgAQHh4uciVERETUU1qtFkql0mwbiWBJBHIgBoMBp0+fhpeXFyQSiVVfW6PRIDw8HKWlpVAoFFZ9bVs2UPsNsO8Dse8Dtd8A+z4Q+25L/RYEAVqtFiEhIZBKza+qGXAjN1KpFGFhYX36HgqFQvRfAjEM1H4D7PtA7PtA7TfAvg/EvttKvy81YtOJC4qJiIjIoTDcEBERkUNhuLEimUyGl19+GTKZTOxS+tVA7TfAvg/Evg/UfgPs+0Dsu732e8AtKCYiIiLHxpEbIiIicigMN0RERORQGG6IiIjIoTDcEBERkUNhuLGS9957D4MGDYKbmxvGjx+P3bt3i11Sr7zyyiuQSCQmtyFDhhgfb2lpwcMPPwxfX1/I5XLcfPPNOHPmjMlrlJSU4LrrroOHhwcCAgLw7LPPor29vb+7ckm//fYbZs+ejZCQEEgkEnz11VcmjwuCgMWLFyM4OBju7u5ISUnBiRMnTNrU1tbirrvugkKhgLe3N+677z40NDSYtDl48CCmTJkCNzc3hIeH48033+zrrl3Spfo+f/78i34PrrnmGpM29tj3jIwMjB07Fl5eXggICMDcuXORn59v0sZav+PZ2dkYM2YMZDIZYmNjsWrVqr7uXrcs6ff06dMv+swXLlxo0sbe+g0AmZmZSExMNF6MLjk5GRs3bjQ+7oifd6dL9d0hP3OBem3t2rWCq6ur8NFHHwmHDx8WFixYIHh7ewtnzpwRu7TL9vLLLwsJCQlCRUWF8Xb27Fnj4wsXLhTCw8OFTZs2CXv37hUmTJggTJw40fh4e3u7MHz4cCElJUXIyckRfvjhB8HPz09IT08Xoztm/fDDD8ILL7wgrF+/XgAgbNiwweTxpUuXCkqlUvjqq6+EAwcOCDfccIMQFRUlNDc3G9tcc801wsiRI4WdO3cKW7duFWJjY4U77rjD+LharRYCAwOFu+66Szh06JDw2WefCe7u7sKKFSv6q5tdulTf582bJ1xzzTUmvwe1tbUmbeyx71dffbXwn//8Rzh06JCQm5srXHvttUJERITQ0NBgbGON3/GTJ08KHh4ewlNPPSUcOXJEePfddwUnJyfhxx9/7Nf+drKk39OmTRMWLFhg8pmr1Wrj4/bYb0EQhG+++Ub4/vvvhePHjwv5+fnC888/L7i4uAiHDh0SBMExP+9Ol+q7I37mDDdWMG7cOOHhhx82/qzX64WQkBAhIyNDxKp65+WXXxZGjhzZ5WP19fWCi4uLsG7dOuN9R48eFQAIO3bsEASh40tTKpUKlZWVxjaZmZmCQqEQdDpdn9beG3/8gjcYDEJQUJCwbNky43319fWCTCYTPvvsM0EQBOHIkSMCAGHPnj3GNhs3bhQkEolQXl4uCIIgvP/++4KPj49J35977jkhPj6+j3tkue7CzZw5c7p9jqP0vaqqSgAgbNmyRRAE6/2O/9///Z+QkJBg8l633XabcPXVV/d1lyzyx34LQscX3eOPP97tcxyh3518fHyElStXDpjP+0KdfRcEx/zMOS3VS62trdi3bx9SUlKM90mlUqSkpGDHjh0iVtZ7J06cQEhICKKjo3HXXXehpKQEALBv3z60tbWZ9HnIkCGIiIgw9nnHjh0YMWIEAgMDjW2uvvpqaDQaHD58uH870gtFRUWorKw06atSqcT48eNN+urt7Y0rrrjC2CYlJQVSqRS7du0ytpk6dSpcXV2Nba6++mrk5+ejrq6un3pzebKzsxEQEID4+HgsWrQINTU1xsccpe9qtRoAoFKpAFjvd3zHjh0mr9HZxlb+bfhjvzv997//hZ+fH4YPH4709HQ0NTUZH3OEfuv1eqxduxaNjY1ITk4eMJ83cHHfOznaZz7gDs60turqauj1epMPHQACAwNx7NgxkarqvfHjx2PVqlWIj49HRUUFlixZgilTpuDQoUOorKyEq6srvL29TZ4TGBiIyspKAEBlZWWXfyedj9mLzlq76suFfQ0ICDB53NnZGSqVyqRNVFTURa/R+ZiPj0+f1N9b11xzDW666SZERUWhsLAQzz//PFJTU7Fjxw44OTk5RN8NBgOeeOIJTJo0CcOHDzfWZY3f8e7aaDQaNDc3w93dvS+6ZJGu+g0Ad955JyIjIxESEoKDBw/iueeeQ35+PtavXw/Avvudl5eH5ORktLS0QC6XY8OGDRg2bBhyc3Md/vPuru+AY37mDDfUpdTUVOOfExMTMX78eERGRuJ///ufqP8Dpf51++23G/88YsQIJCYmIiYmBtnZ2ZgxY4aIlVnPww8/jEOHDmHbtm1il9Kvuuv3Aw88YPzziBEjEBwcjBkzZqCwsBAxMTH9XaZVxcfHIzc3F2q1Gl988QXmzZuHLVu2iF1Wv+iu78OGDXPIz5zTUr3k5+cHJyeni1bVnzlzBkFBQSJVZX3e3t4YPHgwCgoKEBQUhNbWVtTX15u0ubDPQUFBXf6ddD5mLzprNff5BgUFoaqqyuTx9vZ21NbWOtzfR3R0NPz8/FBQUADA/vv+yCOP4LvvvsPmzZsRFhZmvN9av+PdtVEoFKL+n4Tu+t2V8ePHA4DJZ26v/XZ1dUVsbCySkpKQkZGBkSNHYvny5Q7/eQPd970rjvCZM9z0kqurK5KSkrBp0ybjfQaDAZs2bTKZz7R3DQ0NKCwsRHBwMJKSkuDi4mLS5/z8fJSUlBj7nJycjLy8PJMvvp9//hkKhcI4FGoPoqKiEBQUZNJXjUaDXbt2mfS1vr4e+/btM7b59ddfYTAYjP9IJCcn47fffkNbW5uxzc8//4z4+HjRp2V6oqysDDU1NQgODgZgv30XBAGPPPIINmzYgF9//fWiaTNr/Y4nJyebvEZnG7H+bbhUv7uSm5sLACafub31uzsGgwE6nc5hP29zOvveFYf4zEVZxuxg1q5dK8hkMmHVqlXCkSNHhAceeEDw9vY2WVlub55++mkhOztbKCoqErZv3y6kpKQIfn5+QlVVlSAIHdsmIyIihF9//VXYu3evkJycLCQnJxuf37l1cNasWUJubq7w448/Cv7+/ja5FVyr1Qo5OTlCTk6OAEB46623hJycHOHUqVOCIHRsBff29ha+/vpr4eDBg8KcOXO63Ao+evRoYdeuXcK2bduEuLg4k+3Q9fX1QmBgoHD33XcLhw4dEtauXSt4eHiIvhXcXN+1Wq3wzDPPCDt27BCKioqEX375RRgzZowQFxcntLS0GF/DHvu+aNEiQalUCtnZ2SbbX5uamoxtrPE73rk99tlnnxWOHj0qvPfee6Juj71UvwsKCoRXX31V2Lt3r1BUVCR8/fXXQnR0tDB16lTja9hjvwVBEP7yl78IW7ZsEYqKioSDBw8Kf/nLXwSJRCL89NNPgiA45ufdyVzfHfUzZ7ixknfffVeIiIgQXF1dhXHjxgk7d+4Uu6Reue2224Tg4GDB1dVVCA0NFW677TahoKDA+Hhzc7Pw0EMPCT4+PoKHh4dw4403ChUVFSavUVxcLKSmpgru7u6Cn5+f8PTTTwttbW393ZVL2rx5swDgotu8efMEQejYDv7SSy8JgYGBgkwmE2bMmCHk5+ebvEZNTY1wxx13CHK5XFAoFMI999wjaLVakzYHDhwQJk+eLMhkMiE0NFRYunRpf3WxW+b63tTUJMyaNUvw9/cXXFxchMjISGHBggUXhXZ77HtXfQYg/Oc//zG2sdbv+ObNm4VRo0YJrq6uQnR0tMl79LdL9bukpESYOnWqoFKpBJlMJsTGxgrPPvusyTVPBMH++i0IgnDvvfcKkZGRgqurq+Dv7y/MmDHDGGwEwTE/707m+u6on7lEEASh/8aJiIiIiPoW19wQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNwQkejmz5+PuXPn9vv7rlq1ChKJBBKJBE888YRFz5k/f77xOV999VWf1kdEl8dZ7AKIyLFJJBKzj7/88stYvnw5xLpYukKhQH5+Pjw9PS1qv3z5cixdutR4qCAR2R6GGyLqUxUVFcY/f/7551i8eDHy8/ON98nlcsjlcjFKA9ARvoKCgixur1QqoVQq+7AiIuotTksRUZ8KCgoy3pRKpTFMdN7kcvlF01LTp0/Ho48+iieeeAI+Pj4IDAzEBx98gMbGRtxzzz3w8vJCbGwsNm7caPJehw4dQmpqKuRyOQIDA3H33Xejurq6xzW///77iIuLg5ubGwIDA3HLLbf09q+BiPoRww0R2aTVq1fDz88Pu3fvxqOPPopFixbh1ltvxcSJE7F//37MmjULd999N5qamgAA9fX1uOqqqzB69Gjs3bsXP/74I86cOYM//elPPXrfvXv34rHHHsOrr76K/Px8/Pjjj5g6dWpfdJGI+ginpYjIJo0cORIvvvgiACA9PR1Lly6Fn58fFixYAABYvHgxMjMzcfDgQUyYMAH/+te/MHr0aLz++uvG1/joo48QHh6O48ePY/DgwRa9b0lJCTw9PXH99dfDy8sLkZGRGD16tPU7SER9hiM3RGSTEhMTjX92cnKCr68vRowYYbwvMDAQAFBVVQUAOHDgADZv3mxcwyOXyzFkyBAAQGFhocXvO3PmTERGRiI6Ohp33303/vvf/xpHh4jIPjDcEJFNcnFxMflZIpGY3Ne5C8tgMAAAGhoaMHv2bOTm5prcTpw40aNpJS8vL+zfvx+fffYZgoODsXjxYowcORL19fW97xQR9QtOSxGRQxgzZgy+/PJLDBo0CM7OvfunzdnZGSkpKUhJScHLL78Mb29v/Prrr7jpppusVC0R9SWO3BCRQ3j44YdRW1uLO+64A3v27EFhYSGysrJwzz33QK/XW/w63333Hd555x3k5ubi1KlT+Pjjj2EwGBAfH9+H1RORNTHcEJFDCAkJwfbt26HX6zFr1iyMGDECTzzxBLy9vSGVWv5Pnbe3N9avX4+rrroKQ4cOxb///W989tlnSEhI6MPqiciaJIJYlwUlIhLZqlWr8MQTT1zWehqJRIINGzaIcmwEEZnHkRsiGtDUajXkcjmee+45i9ovXLhQ1CsqE9GlceSGiAYsrVaLM2fOAOiYjvLz87vkc6qqqqDRaAAAwcHBFp9JRUT9h+GGiIiIHAqnpYiIiMihMNwQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FD+P1bdBYJlA0RMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load model\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# Step 1: Parameter replacement\n",
    "parameter_values = model.default_parameter_values\n",
    "parameter_values.process_model(model)\n",
    "geometry = model.default_geometry\n",
    "parameter_values.process_geometry(geometry)\n",
    "\n",
    "# Step 2: Discretisation\n",
    "mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)\n",
    "disc = pybamm.Discretisation(mesh, model.default_spatial_methods)\n",
    "disc.process_model(model)\n",
    "\n",
    "# Step 3: Solver setup\n",
    "#    The `set_up` method is normally called by the `solve` method if not called explicitly.\n",
    "#    We also need to update `_model_set_up` to include the initial conditions or the solver\n",
    "#    will call `set_up` again during `solve`.\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "solver.set_up(model)\n",
    "solver._model_set_up.update(\n",
    "    {model: {\"initial conditions\": model.concatenated_initial_conditions}}\n",
    ")\n",
    "\n",
    "# Step 4: Solver solve\n",
    "solution = solver.solve(model, [0, 3600])\n",
    "\n",
    "# Step 5: Post-processing\n",
    "t_evals = np.linspace(0, 3600, 100)\n",
    "voltage = solution[\"Terminal voltage [V]\"](t_evals)\n",
    "\n",
    "# Plot the results\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(t_evals, voltage)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Voltage [V]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll instrument this script to time each part of the pipeline, using the high-resolution timer provided by the `time` module in Python. This will allow us to see how long each part of the pipeline takes, and where the performance bottlenecks are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter replacement took 0.073 seconds\n",
      "Discretisation took 0.238 seconds\n",
      "Solver setup took 0.076 seconds\n",
      "Solver solve took 0.028 seconds\n",
      "Post-processing took 0.002 seconds\n",
      "Total time taken: 0.418 seconds\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "# load model\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# Step 1: Parameter replacement\n",
    "total_start = time.perf_counter()\n",
    "start = time.perf_counter()\n",
    "parameter_values = model.default_parameter_values\n",
    "parameter_values.process_model(model)\n",
    "geometry = model.default_geometry\n",
    "parameter_values.process_geometry(geometry)\n",
    "end = time.perf_counter()\n",
    "print(f\"Parameter replacement took {end - start:.3f} seconds\")\n",
    "\n",
    "# Step 2: Discretisation\n",
    "start = time.perf_counter()\n",
    "mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)\n",
    "disc = pybamm.Discretisation(mesh, model.default_spatial_methods)\n",
    "disc.process_model(model)\n",
    "end = time.perf_counter()\n",
    "print(f\"Discretisation took {end - start:.3f} seconds\")\n",
    "\n",
    "# Step 3: Solver setup\n",
    "start = time.perf_counter()\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "solver.set_up(model)\n",
    "solver._model_set_up.update(\n",
    "    {model: {\"initial conditions\": model.concatenated_initial_conditions}}\n",
    ")\n",
    "end = time.perf_counter()\n",
    "print(f\"Solver setup took {end - start:.3f} seconds\")\n",
    "\n",
    "# Step 4: Solver solve\n",
    "start = time.perf_counter()\n",
    "solution = solver.solve(model, [0, 3600])\n",
    "end = time.perf_counter()\n",
    "print(f\"Solver solve took {end - start:.3f} seconds\")\n",
    "\n",
    "# Step 5: Post-processing\n",
    "start = time.perf_counter()\n",
    "t_evals = np.linspace(0, 3600, 100)\n",
    "voltage = solution[\"Terminal voltage [V]\"](t_evals)\n",
    "end = time.perf_counter()\n",
    "print(f\"Post-processing took {end - start:.3f} seconds\")\n",
    "total_end = time.perf_counter()\n",
    "print(f\"Total time taken: {total_end - total_start:.3f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this output, we can see that the actual time taken to solve the equations is only a small percentage (less than 10%) of the total time taken to run the simulation. In fact, the most time-consuming step is the discretisation, which takes roughly 50% of the total time to discretise the spatial derivatives of the model. Even when this has all been done, the solver setup takes more than double the time of the actual solve. All of these numbers emphasises the fact that there is significant overhead in setting up the equations and the solver before the actual solution can be computed. Therefore, if you are looking to optimise the performance of your simulations, you should take care to ensure that you minimise the time spent in the setup steps of the pipeline, particularly if you are running many simulations with the same model and/or discretisation.\n",
    "\n",
    "Note you can also use the logging features in PyBaMM to automatically print out the time taken for the solver setup and solve time. This can be done by setting the logging level to `INFO` before running the simulation, like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-01-07 17:32:33.236 - [INFO] base_model._build_model(829): Start building Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.268 - [INFO] base_battery_model.build_model(1072): Finish building Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.269 - [INFO] parameter_values.process_model(463): Start setting parameters for Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.340 - [INFO] parameter_values.process_model(532): Finish setting parameters for Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.344 - [INFO] discretisation.process_model(142): Start discretising Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.581 - [INFO] discretisation.process_model(244): Finish discretising Doyle-Fuller-Newman model\n",
      "2025-01-07 17:32:33.582 - [INFO] base_solver.solve(758): Start solving Doyle-Fuller-Newman model with IDA KLU solver\n",
      "2025-01-07 17:32:33.584 - [INFO] base_solver.set_up(138): Start solver set-up\n",
      "2025-01-07 17:32:33.637 - [INFO] base_solver.set_up(280): Finish solver set-up\n",
      "2025-01-07 17:32:33.689 - [INFO] base_solver.solve(979): Finish solving Doyle-Fuller-Newman model (the solver successfully reached the end of the integration interval)\n",
      "2025-01-07 17:32:33.690 - [INFO] base_solver.solve(980): Set-up time: 83.762 ms, Solve time: 21.137 ms (of which integration time: 20.928 ms), Total time: 104.899 ms\n"
     ]
    }
   ],
   "source": [
    "pybamm.set_logging_level(\"INFO\")\n",
    "\n",
    "# load model\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# Step 1: Parameter replacement\n",
    "parameter_values = model.default_parameter_values\n",
    "parameter_values.process_model(model)\n",
    "geometry = model.default_geometry\n",
    "parameter_values.process_geometry(geometry)\n",
    "\n",
    "# Step 2: Discretisation\n",
    "var = pybamm.standard_spatial_vars\n",
    "mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)\n",
    "disc = pybamm.Discretisation(mesh, model.default_spatial_methods)\n",
    "disc.process_model(model)\n",
    "\n",
    "# Step 3 & 4: Solver setup\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "t_eval = [0, 3600]\n",
    "solution = solver.solve(model, t_eval)\n",
    "\n",
    "# Step 5: Post-processing\n",
    "t_interp = np.linspace(0, 3600, 100)\n",
    "voltage = solution[\"Terminal voltage [V]\"](t_interp)\n",
    "\n",
    "pybamm.set_logging_level(\"WARNING\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Simulation Class\n",
    "\n",
    "The `pybamm.Simulation` class is a convenience class that wraps up the entire pipeline into a single class that can be used to setup the model equations and solver and run the simulation in a single call. This class is an easy-to-use wrapper for the PyBaMM pipeline, and provides many useful features for initialising simulations and running experiments. However, it is important to understand that the `pybamm.Simulation` class still performs all the steps of the pipeline, and so the performance characteristics of the `pybamm.Simulation` class are the same as the individual steps of the pipeline.\n",
    "\n",
    "Lets see how we can use the `pybamm.Simulation` class to run the same simulation as above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLEklEQVR4nO3dd3zT1f4/8FfSkY40adO96aJAoYzKKFspaFUE13VdCw4U3PPrrQPFe7Uo93pFr/byE72gXsSLghurIkVANi2UVWhp6aCldCXpStvk8/ujNBBpQ0rTfpL09Xw88pAmJ8n7kF7yumd8jkQQBAFEREREDkIqdgFERERE1sRwQ0RERA6F4YaIiIgcCsMNERERORSGGyIiInIoDDdERETkUBhuiIiIyKEw3BAREZFDYbghIiIih8JwQ0RERA5F1HCTmZmJxMREKBQKKBQKJCcnY+PGjWaf8/bbbyM+Ph7u7u4IDw/Hk08+iZaWln6qmIiIiGyds5hvHhYWhqVLlyIuLg6CIGD16tWYM2cOcnJykJCQcFH7NWvW4C9/+Qs++ugjTJw4EcePH8f8+fMhkUjw1ltvidADIiIisjUSWzs4U6VSYdmyZbjvvvsueuyRRx7B0aNHsWnTJuN9Tz/9NHbt2oVt27ZZ9PoGgwGnT5+Gl5cXJBKJ1eomIiKiviMIArRaLUJCQiCVmp94EnXk5kJ6vR7r1q1DY2MjkpOTu2wzceJEfPrpp9i9ezfGjRuHkydP4ocffsDdd9/d7evqdDrodDrjz+Xl5Rg2bJjV6yciIqK+V1pairCwMLNtRA83eXl5SE5ORktLC+RyOTZs2NBt+LjzzjtRXV2NyZMnQxAEtLe3Y+HChXj++ee7ff2MjAwsWbLkovtLS0uhUCis1g8iIiLqOxqNBuHh4fDy8rpkW9GnpVpbW1FSUgK1Wo0vvvgCK1euxJYtW7oMONnZ2bj99tvxt7/9DePHj0dBQQEef/xxLFiwAC+99FKXr//HkZvOvxy1Ws1wQ0REZCc0Gg2USqVF39+ih5s/SklJQUxMDFasWHHRY1OmTMGECROwbNky432ffvopHnjgATQ0NFxyDg7o2V8OERER2YaefH/b3HVuDAaDyUjLhZqami4KME5OTgA6FhoRERERibrmJj09HampqYiIiIBWq8WaNWuQnZ2NrKwsAEBaWhpCQ0ORkZEBAJg9ezbeeustjB492jgt9dJLL2H27NnGkENEREQDm6jhpqqqCmlpaaioqIBSqURiYiKysrIwc+ZMAEBJSYnJSM2LL74IiUSCF198EeXl5fD398fs2bPx2muvidUFIiIisjE2t+amr3HNDRERkf2x6zU3RERERL3BcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNxY0en6Zhyr1IhdBhER0YDGcGMlG/MqMG3ZZqSvz+M5V0RERCJiuLGSpEE+kEokyCmpx/aCGrHLISIiGrAYbqwkwMsNd46PAAAs33ScozdEREQiYbixooXTYuDqLMWe4jrsPFkrdjlEREQDEsONFQUq3HDbFeEAgHc2nRC5GiIiooGJ4cbKFk6PgYuTBDtO1mBPMUdviIiI+hvDjZWFervjliSO3hAREYmF4aYPPDQ9Bs5SCbaeqMb+kjqxyyEiIhpQGG76QLjKAzeNCQUAvMvRGyIion7FcNNHHpoeC6kE2Jx/FluOnxW7HCIiogGD4aaPDPLzxC1JYQCA+1fvwfr9ZSJXRERENDAw3PShV+cMx3UjgtGmF/DU/w7g7V94cT8iIqK+xnDTh9xcnPDuHaOxcFoMAODtX07g6XUH0NpuELkyIiIix8Vw08ekUgn+kjoEr984Ak5SCdbvL8cdH+xE4dkGsUsjIiJySAw3/eTO8RH4aP5YyGXO2HeqDqnLt+LdTSc4ikNERGRlDDf9aNpgf2x8fAqmDvZHa7sB//j5OGa/u43XwiEiIrIihpt+Fq7ywOp7xuLt20ZB5emK/DNa3Jz5O5774iCqNC1il0dERGT3GG5EIJFIMHd0KDY9NQ03jwmDIACf7y3F9L9nY/kvJ9DU2i52iURERHZLIgywvckajQZKpRJqtRoKhULscgAA+07V4m/fH0VOST0AIFAhwzOz4nHTmDA4SSXiFkdERGQDevL9zXBjIwRBwHcHK/DGj8dQVtcMAIgNkOOZWYNxdUIQJBKGHCIiGrgYbsyw1XDTqaVNj9W/F+P97EKom9sAAIlhSjx7dTwmx/ox5BAR0YDEcGOGrYebTpqWNnzw20l8uK0ITa16AMC4KBWemjkYE6J9Ra6OiIiofzHcmGEv4aZTdYMO728uxKc7T6FV33FNnORoXzw5czDGRalEro6IiKh/MNyYYW/hplOFuhnvby7E2j0laNN3fGSTY/3w6FWxGM+RHCIicnAMN2bYa7jpVF7fjPc2F2Dd3lJjyBkfpcJjM+IwMcaXa3KIiMghMdyYYe/hplNZXRMyswuxbm+ZcbpqTIQ3Hr0qDtPj/RlyiIjIoTDcmOEo4aZThboZK7acxGe7S6A7d07V0GAFHr4yBqnDg3mdHCIicggMN2Y4WrjpVKVpwQdbT+K/u0qMu6ui/DyxcFo05o4OhczZSeQKiYiILh/DjRmOGm461TW2YtXvxVj1e7HxOjmBChnunRSFO8dHwMvNReQKiYiIeo7hxgxHDzedGnTtWLPrFD7cVoQzGh0AwEvmjLsmROLeSYMQoHATuUIiIiLLMdyYMVDCTSddux5f557Gii2FKDzbCABwcZLghpGhWDA1CkOCHP/vgIiI7B/DjRkDLdx0MhgEbDpWhX9vKcS+U3XG+6fE+eH+KdGYGsejHYiIyHYx3JgxUMPNhfaX1GHl1pP48VAlDOc+/bgAOe6dHIW5o0Lh7srFx0REZFsYbsxguDmvtLYJH24rwrq9pWg8t8PK28MFd46LwN3JkQhWuotcIRERUQeGGzMYbi6maWnD//aUYtXvxSirawYAOEkluCYhCGnJkRgXpeKUFRERiYrhxgyGm+7pDQJ+PnIGH20vwu6iWuP9Q4K8MG/iIMwZFQIPV2cRKyQiooGK4cYMhhvLHK3Q4OMdxdiQU46Wto4rH3u5OePmMWH484RIxAbIRa6QiIgGkp58f0v7qaYuZWZmIjExEQqFAgqFAsnJydi4cWO37adPnw6JRHLR7brrruvHqgeGocEKZNyUiF3pKXjxuqGIUHlA29KOVb8XI+WtLbjj/+3E9wcr0HruyAciIiJbIerIzbfffgsnJyfExcVBEASsXr0ay5YtQ05ODhISEi5qX1tbi9bWVuPPNTU1GDlyJFauXIn58+db9J4cubk8BoOArQXV+HTnKWw6esa4y8pP7opbksJx+9hwDPLzFLdIIiJyWHY9LaVSqbBs2TLcd999l2z79ttvY/HixaioqICnp2VfrAw3vVde34y1u0vw+Z5SVGl1xvsnxfri9rERmJUQyLOsiIjIqnry/W0zq0P1ej3WrVuHxsZGJCcnW/ScDz/8ELfffrvZYKPT6aDTnf8C1mg0va51oAv1dsfTs+Lx2Iw4/HqsCp/tLsGW42exvaAG2wtq4OPhghtHh+G2seGID/ISu1wiIhpgRB+5ycvLQ3JyMlpaWiCXy7FmzRpce+21l3ze7t27MX78eOzatQvjxo3rtt0rr7yCJUuWXHQ/R26sq6yuCZ/vKcW6vWWo1LQY7x8Z7o0/XRGG2SNDoOChnUREdJnsalqqtbUVJSUlUKvV+OKLL7By5Ups2bIFw4YNM/u8Bx98EDt27MDBgwfNtutq5CY8PJzhpo/oDQJ+O34Wn+8pxS9Hz6D93OIcmbMU1wwPwq1J4ZgY4wuplNfNISIiy9lVuPmjlJQUxMTEYMWKFd22aWxsREhICF599VU8/vjjPXp9rrnpP2e1OnyVU451+0px/EyD8f4QpRvmjg7FzUlhiPHnlnIiIro0u1xz08lgMJiMtHRl3bp10Ol0+POf/9xPVdHl8PeSYcHUaNw/JQoHy9T4Yl8Zvs4tx2l1C97PLsT72YUYGe6NW8aE4vrEEPh4uopdMhEROQBRR27S09ORmpqKiIgIaLVarFmzBm+88QaysrIwc+ZMpKWlITQ0FBkZGSbPmzJlCkJDQ7F27doevydHbsTV0qbHpqNVWL+/DNnHz0J/btrKxUmC6fEBuHlMKK4cEsDdVkREZMJuRm6qqqqQlpaGiooKKJVKJCYmGoMNAJSUlEAqNb3OYH5+PrZt24affvpJjJKpl9xcnHBdYjCuSwzGWa0OX+eWY0NOOQ6f1uDnI2fw85EzULq74NoRwZg7KgRjB6m4PoeIiHrE5tbc9DWO3Nim/Eot1ueU4euc0ya7rUKUbpg9KgRzR4ViSJAXD/AkIhqg7HpBcV9juLFteoOAnSdr8HVuOTbmVUKrazc+Fhcgx+yRIbg+MRjRXIhMRDSgMNyYwXBjP1ra9MjOr8JXOafxa36VyTlWw0MVuD4xBNeNCEa4ykPEKomIqD8w3JjBcGOfNC1t+OnwGXx74DS2FVQbFyIDwKhwb1yfGIxrRwQjxNtdxCqJiKivMNyYwXBj/2oadNh4qBLfHTyNXUW1uPA3eHSEN64bEYxrhgchzIcjOkREjoLhxgyGG8dSpW3Bj4cq8d3BCuwpNg06I8OUSB0RjGsSgnhiORGRnWO4MYPhxnGd0bQg63AlfsirwO6iWlwwc4UhQV6YlRCEaxKCMDSYu66IiOwNw40ZDDcDw1mtDlmHK5F1uBI7CmuMZ1wBQLjKHbOGBWHWsEAkRfrA2Ulq5pWIiMgWMNyYwXAz8NQ3tWLT0Sr8eLgSvx0/C90Fu658PFwwY2ggUoYGYkqcHzxlNnciCRERgeHGLIabga2ptR2/Ha/GT0cqseloFdTNbcbHXJ2lmBjjixlDAzFjSAB3XhER2RCGGzMYbqhTu96A3cW1+OVIFX45egYltU0mjw8NVuCqIf64akggRoV7w4nHQBARiYbhxgyGG+qKIAgoqGrAz0fPYNPRKuwvqTPZeaXydMXUOD9cOSQAU+L8oeIJ5kRE/YrhxgyGG7JEbWMrthyvwqajVdhy/Cy0LeePgZBIgJFh3rgyPgDT4v0xIlTJUR0ioj7GcGMGww31VJvegP2n6pB9/Cw2H6vCsUqtyeM+Hi6YEuePaYP9MSXODwEKN5EqJSJyXAw3ZjDcUG9VqJuxJf8sthw/i20nqk0O9wQ6rqkz9VzQGTtIBTcXJ5EqJSJyHAw3ZjDckDW16Q3ILa03hp1Dp9Uma3VkzlKMHaTCpFg/TI71w7AQBaewiIguA8ONGQw31JdqGnTYXliDrcfPYuuJalRqWkwe9/ZwQXK0LybG+mFijC+i/Tx5tWQiIgsw3JjBcEP9RRAEFJ5twLYT1dhWUIOdJ2vQ8IcprGClG5JjfI2BJ5TX1iEi6hLDjRkMNySWdr0BB8rU+L2gGtsLq7H/VD1a9QaTNhEqD0yM8UVyjC/GR/kiSMnFyUREAMONWQw3ZCuaW/XYe6oWOwpr8HthDfLK1dAbTP/nGOXniQnRKkyIZtghooGN4cYMhhuyVdqWNuwp7gg7u4pqcahcjT9kHUT6emB8lArjo3wxPlqFMB8PcYolIupnDDdmMNyQvVA3t2FvcS12nuw+7IR6u2NclArjo1QYF6VCFBcoE5GDYrgxg+GG7JWmpQ37iuuws6gGu07WdjmN5Sd3xdhBHUFn7CAVhgZz6zkROQaGGzMYbshRNOrakVNSj91FNdhZVIvc0nq0tpsuUPaSOSNpkA/GRakwbpAKI8KUkDnzooJEZH8YbsxguCFHpWvX42CZGruLarG7qBb7TtVdtPVc5izFqHBvjI9SYWyUCmMifOApcxapYiIiyzHcmMFwQwNFu96AY5Va7CqqxZ6iWuwurkVtY6tJGyepBMNDlecWKatwxSAVlO4uIlVMRNQ9hhszGG5ooOq8qODuojrsKe4Y3SmvbzZpI5EAQ4MUmBDtiwnRHWt3vD1cRaqYiOg8hhszGG6Iziura8Ke4lrsOtkRdk5WN5o83hl2kmN8MSnWF+OifCHnNBYRiYDhxgyGG6LuVWlasKuoY/v5zpM1KDxrGnacpRKMDPfGpBhfTB3sj1Hh3nB2kopULRENJAw3ZjDcEFmuStuCXSdr8XthDbYXVKOktsnkcYWbM6YM9se0wf6YHu+PAC9eQZmI+gbDjRkMN0SXr7S2CdsLqrG1oBrbTlRD3dxm8vjoCG9cnRCEqxOCEOXnKVKVROSIGG7MYLghso7Og0C35Fch+/hZHCxTmzweFyDHtSOCccOoEMT4y0WqkogcBcONGQw3RH2jUt2Cn4+ewU+HK7GjsAbtF1w9eXioAjeMDMHskSEIVrqLWCUR2SuGGzMYboj6nrq5DZuOnsG3B07jtxPVxmMiJBJgcqwfbkkKw9UJQXBz4dWSicgyDDdmMNwQ9a+aBh1+OFSJb3LLsae4zni/l5szZo8Mwa1JYRgV7s0DP4nILIYbMxhuiMRzqqYRX+4vx5f7ykwuIBgXIMdtY8Mxd3Qo/OQyESskIlvFcGMGww2R+AwGATtP1mDdvjL8kFcB3bkDP52lEswYGoBbksIxPd4fLryGDhGdw3BjBsMNkW3RtLTh2wOn8b89pThwwY4rX09XzBkViluSwjAshP9bJRroGG7MYLghsl3HKjX4cl8ZNuScRnWDznj/kCAvzB4ZgusTgxHpy+vnEA1EDDdmMNwQ2b52vQG/nTiLL/aV4ZcjVWjVG4yPjQxT4vrEEFwzPAjhKg8RqySi/sRwYwbDDZF9UTe1IetwJb49eBrbC6pxweVzMCTICylDAzFzWCBGhCohlXLHFZGjYrgxg+GGyH6d1erw46EKfJ9XgT3Fdcbr5wCAv5cMk2P9MCXOD5Nj/RCg4DlXRI6E4cYMhhsix1Df1IrN+VX4+cgZbMk/i8ZWvcnj8YFemBCtwtgoFcYNUjHsENk5hhszGG6IHI+uXY99p+qw7UQ1tp6oxqHTavzxX7ZBvh5IilRhVIQ3Rod7Iz7Ii1vNiewIw40ZDDdEjq+2sRU7Cmuwu6gGu4vrcKxSc1HYkTlLMTxUiRGhSgwLUSAhRIG4AC+4OjPwENkihhszGG6IBh51cxv2n6rD/pI65JbW40BpPTQt7Re1c3WSIi5QjvggLwwJ8kJ8kAJDgrwQ4CXj8RBEIrObcJOZmYnMzEwUFxcDABISErB48WKkpqZ2+5z6+nq88MILWL9+PWpraxEZGYm3334b1157rUXvyXBDRAaDgKKaRhworcfh0xocPq3G4dMaaLsIPACgdHdBbIAccQFyxJ67xfjLEeLtDifu0CLqF3YTbr799ls4OTkhLi4OgiBg9erVWLZsGXJycpCQkHBR+9bWVkyaNAkBAQF4/vnnERoailOnTsHb2xsjR4606D0ZboioK4IgoLS2GUcqNMiv1CL/jAbHKrUorm402X5+IVdnKaJ8PRHt74koP08M8vPEIF9PDPLzgL+coz1E1mQ34aYrKpUKy5Ytw3333XfRY//+97+xbNkyHDt2DC4uLpf1+gw3RNQTLW16FFU34kRVAwrOaDv+W9WAUzVNJhcX/CNPVyeEqzwQofJApG/Hf8NUHgj38UCYjzvcXJz6sRdE9s8uw41er8e6deswb9485OTkYNiwYRe1ufbaa6FSqeDh4YGvv/4a/v7+uPPOO/Hcc8/Byanrfyh0Oh10uvOXcddoNAgPD2e4IaJe0RsElNc1o7C6ASfPNqK4uhHFNR238rrmbkd7OvnJZQhXuSPE2x1h3u4I9XFHqHfHzyFKdyjcnTnyQ3SBnoQb536qqVt5eXlITk5GS0sL5HI5NmzY0GWwAYCTJ0/i119/xV133YUffvgBBQUFeOihh9DW1oaXX365y+dkZGRgyZIlfdkFIhqAnKQSRPh6IMLXA1fGmz6ma9ejrK4ZJbVNKKlpQkltE07VNKGsrgmltU1obNWjukGH6gYdckrqu3x9D1cnBCvdEOLtjmClG4KUnf91Q7DSDcEKBiCi7og+ctPa2oqSkhKo1Wp88cUXWLlyJbZs2dJlwBk8eDBaWlpQVFRkHKl56623sGzZMlRUVHT5+hy5ISJbIggC1M1tKK1tRmldE07XN6O8vhnldR3/rVC3oLax1aLXcndxQpDSDUEKN5Pg0xmEgpVuUHm6MgCRQ7CrkRtXV1fExsYCAJKSkrBnzx4sX74cK1asuKhtcHAwXFxcTKaghg4disrKSrS2tsLV1fWi58hkMshksr7rABFRD0gkEnh7uMLbwxUjwpRdtmlp06NC3YLT9c04Xd+MSnULKjQtHf9Vt6BC3Yz6pjY0n1sPVFTd2O37uTpLjUEnxPuCqS9vd4R6d9zn4Sr6VwGRVdncb7TBYDAZabnQpEmTsGbNGhgMBkilHRfaOn78OIKDg7sMNkRE9sjNxQlRfh07sLrT0qbHGU1H2KlUt6DyXPg5Xd+MSk0LTte3oLpBh9Z2A07VdEyLdUfl6YrQc8EnzKfjFnpu4XOYjzu83C5vAweRWEQNN+np6UhNTUVERAS0Wi3WrFmD7OxsZGVlAQDS0tIQGhqKjIwMAMCiRYvwr3/9C48//jgeffRRnDhxAq+//joee+wxMbtBRNTv3FycEOnriUjf7gNQa7sBZzQdgadC3YLT6uZzo0EtxukwbUs7ahtbUdvYirxydZevo3R3QbjKHeE+HghXeRh3gUWoPBDq7c6rOpPNETXcVFVVIS0tDRUVFVAqlUhMTERWVhZmzpwJACgpKTGO0ABAeHg4srKy8OSTTyIxMRGhoaF4/PHH8dxzz4nVBSIim+XqLDWGke5oWtpQXteMsrpmlNc1oby+488dtybUNbVB3dwGdXkbDpVrLnq+RAKEKN0RrnJHpMqzY5H1ue3vkb6eULpz1If6n+gLivsbr3NDRGS5Bl37uV1ezSitbULpBX8uqW1Cc5ve7PN9PFzOX9zw3AUOo/3kGOTnweku6hG7vM5Nf2G4ISKyDkEQUN3Q2rHlvbYRp2rOb30vrmlCdUPX6yc7+clliD63tijq3FWeY/zliPT14IntdBGGGzMYboiI+keDrh2najpCT1F1I07VNKK4ugknqxvNBh8nqQSRKg9E+8sRE+CJWP/zZ3pxtGfgYrgxg+GGiEh82pa2c0Gnwbid/eTZRpw824DG1u6nugIVMgwO9MLgQC/EB3ohPsgLcYFybmcfABhuzGC4ISKyXYIg4IxGh8KzDSg823GOV+etStv1aI9EAkSqPDAkSIEhwV4YEqTA0GAvRKg8eAFDB8JwYwbDDRGRfVI3t50LOlocq9Ti+Bkt8isbup3i8nJzxrBgBRJClEgIUSAxTIkYfzmkUgYee8RwYwbDDRGRY6lp0CG/UoujlVocq9DgaKUGxysbujy1XS5zxvBQBUaGeWNUuDeSIn0QoHAToWrqKYYbMxhuiIgcX5vegIKqBhw+rcGhcjUOn1bjULmmy63rYT7uuCLSB0mRPhgX5YvBgXJOZ9kghhszGG6IiAamdr0BBWcbcLBUjdyyeuSU1CO/UgPDH74FfT1dMSHaFxOiVUiO8UOMvyfDjg1guDGD4YaIiDppW9qQW1qPfafqsLe4DvtO1V00uhPq7Y5p8f6YPtgfE2P9IJdxZ5YYGG7MYLghIqLutLYbcLCsHjsKa7DjZA32nqpDa/v5tTsuThJMiPbF1QlBuDohCP5eMhGrHVgYbsxguCEiIks1tbZj18laZOdXIfv4WZPT1SUSYGykCqkjgnB9YgiDTh9juDGD4YaIiC7XybMNyDp8Bj8eqsCBsvOnqDtJJZg+2B83J4VhxtAAyJydRKzSMTHcmMFwQ0RE1lBe34wfD1Xi2wOnkVtab7xf6e6CG0eHYt7EQYjy8xSvQAfDcGMGww0REVlbQVUD1u8vw/r95ajUtADomLa6Kj4A906OwsQYX+646iWGGzMYboiIqK/oDQK2FVTj49+LselYlfH++EAvPHRlDGYnhvAKyZeJ4cYMhhsiIuoPJ882YPXvxVi3rwxN5w4DjQ/0wlOzBmPWsECO5PQQw40ZDDdERNSf1M1t+Pj3Yvy/rSehbWkHACSGKfHMrHhMHewvcnX2g+HGDIYbIiISg7qpDf9vayH+s73YOJJz7YggLL4+AUFKnm91KQw3ZjDcEBGRmKobdHhvcwE+3nEKeoMAucwZT88ajLTkQXDiepxuMdyYwXBDRES24MhpDV74Kg85JfUAgOGhCrxxcyISQpTiFmajevL9Le2nmoiIiOgCw0IU+HLhRLx+4wgo3JxxqFyDG9//HR/vKMYAG3ewOoYbIiIikUilEtw5PgKbnp6OlKGBaG03YPHXh7Ho0/1QN7eJXZ7dYrghIiISmb+XDB+kJWHx9cPg4iTBj4crcd07W5FTUid2aXaJ4YaIiMgGSCQS3Ds5Cl8umogIlQfK6prxpxU78HVuudil2R2GGyIiIhuSGOaN7x6bjGsSgtCmF/D42lys3HpS7LLsCsMNERGRjVG4ueD9u8bg3klRAIC/fX8Ur31/BAYDFxpbguGGiIjIBkmlErx0/VCkpw4BAHywtQhP/S8Xre0GkSuzfQw3RERENkoikeDBaTF4608j4SyV4Kvc03jov/vQrmfAMYfhhoiIyMbdNCYMK+ddAZmzFL8crUL6+jxeC8cMZ0saqVSqHr2oRCLB/v37ERkZeVlFERERkanp8QF4784xePDTfVi3rwz+XjL83zVDxC7LJlkUburr6/H2229Dqbz0JaEFQcBDDz0EvV7f6+KIiIjovJRhgXj9xuF47ss8vJ9diAAvGeafW3RM51kUbgDg9ttvR0BAgEVtH3300csuiIiIiLp329gInNXq8PefjmPJd0fg5yXD9YkhYpdlUywKNwZDzxYuabXayyqGiIiILu3hK2NRpdXh4x2n8OTnuQhWuiMp0kfssmyGxQuKv/vuux6HHCIiIrI+iUSCl2cnGC/098TnOdC28CyqThaHm7lz5yI8PBwvvPACCgoK+rImIiIiugQnqQRv3pqIMB93lNY2Y/HXh8UuyWZYHG6Kiorw4IMPYu3atYiPj8e0adPwySefoLm5uS/rIyIiom4o3Fzw9m2jIJUAG3LKeQ7VORaHm/DwcCxevBiFhYX45ZdfMGjQICxatAjBwcFYuHAh9uzZ05d1EhERUReuGKTCo1fFAQBe3HAIpbVNIlckvsu6iN+VV16J1atXo6KiAsuWLUNeXh4mTJiAkSNHWrs+IiIiuoRHr4pFUqQPtLp2PPF57oC/gnGvrlDs5eWFGTNm4Morr4S3tzeOHDlirbqIiIjIQs5OUrx92yh4yZyx71Qd3ttcKHZJorqscNPc3IyPP/4Y06dPR1xcHNauXYunnnoKxcXFVi6PiIiILBGu8sBf5w4HALy3uWBAT0/1KNzs3LkTDzzwgHGdTVhYGH755RcUFBTghRdeQGhoaF/VSURERJcwZ1QIJsX6olVvwN9/yhe7HNFYHG6GDRuGSZMmYf/+/cjIyEBFRQU+/fRTXHnllX1ZHxEREVlIIpEgPXUoAODr3NM4VK4WuSJxWBxuUlJSsH//fuzduxeLFi2y6JwpIiIi6l/DQ5WYM6rjOIaMjUcH5OnhFoebd955h7uhiIiI7MAzs+Lh6iTF9oIa/HaiWuxy+p1F4WbMmDGoq6uz+EUnT56M8nJeSIiIiEgM4SoPpCVHAgAyfjgKvWFgjd5YdHBmbm4uDhw4AJVKZdGL5ubmQqfTXbJdZmYmMjMzjbusEhISsHjxYqSmpnbZftWqVbjnnntM7pPJZGhpabGoLiIiooHi4Stj8fneUhyr1OKrnHLcnBQmdkn9xqJwAwAzZsyweN5OIpFY1C4sLAxLly5FXFwcBEHA6tWrMWfOHOTk5CAhIaHL5ygUCuTnn18Bbul7ERERDSQ+nq54+MpYLN14DP/4KR/XJQbDzcVJ7LL6hUXhpqioqMcvHBZ26YQ4e/Zsk59fe+01ZGZmYufOnd2GG4lEgqCgoB7XQ0RENNDMnzgIq38vxml1Cz7fU4p5EweJXVK/sCjcREZG9nUd0Ov1WLduHRobG5GcnNxtu4aGBkRGRsJgMGDMmDF4/fXXuw1CAKDT6UymyDQajVXrJiIislVuLk5YND0Gi78+jP/uOoW05MgBMePRq+MXrCEvLw9yuRwymQwLFy7Ehg0bMGzYsC7bxsfH46OPPsLXX3+NTz/9FAaDARMnTkRZWVm3r5+RkQGlUmm8hYeH91VXiIiIbM7c0aFwc5Hi+JkG7C+xfHOQPZMIIm+Ab21tRUlJCdRqNb744gusXLkSW7Zs6TbgXKitrQ1Dhw7FHXfcgb/+9a9dtulq5CY8PBxqtRoKhcJq/SAiIrJVz647gHX7ynDzmDD840/2eVkXjUYDpVJp0fe36CM3rq6uiI2NRVJSEjIyMjBy5EgsX77coue6uLhg9OjRKCgo6LaNTCaDQqEwuREREQ0kd4yPAAB8d/A01E1tIlfT90QPN39kMBgs2kYOdKzTycvLQ3BwcB9XRUREZL9Gh3tjSJAXdO0GbMjpfimHo7iscFNfX4+VK1ciPT0dtbW1AID9+/f3+MJ96enp+O2331BcXIy8vDykp6cjOzsbd911FwAgLS0N6enpxvavvvoqfvrpJ5w8eRL79+/Hn//8Z5w6dQr333//5XSDiIhoQJBIJLjz3OjNZ7tLHf5IBouvc9Pp4MGDSElJgVKpRHFxMRYsWACVSoX169ejpKQEH3/8scWvVVVVhbS0NFRUVECpVCIxMRFZWVmYOXMmAKCkpARS6fn8VVdXhwULFqCyshI+Pj5ISkrC77//btH6HCIiooFszqhQvP7DUeSf0WJ/ST2SIn3ELqnP9HhBcUpKCsaMGYM333wTXl5eOHDgAKKjo/H777/jzjvvNF5t2Fb1ZEESERGRI3lm3QF8sa8MtySF4e+32tfC4j5dULxnzx48+OCDF90fGhqKysrKnr4cERER9ZM7xl2wsLjZcRcW9zjcyGSyLi+Ed/z4cfj7+1ulKCIiIrK+MRHeiA/0QkubAV/lOO4B1z0ONzfccANeffVVtLV1JD6JRIKSkhI899xzuPnmm61eIBEREVmH6cLiEpGr6Ts9Djf/+Mc/0NDQgICAADQ3N2PatGmIjY2Fl5cXXnvttb6okYiIiKxk7uhQuDhJcKxSi6LqRrHL6RM93i2lVCrx888/Y9u2bTh48CAaGhowZswYpKSk9EV9REREZEVKdxdcEanCjpM1yM6vQpRflNglWV2Pw02nyZMnY/LkydashYiIiPrB9Hh/7DhZgy3Hz+KeSQw3eOedd7q8XyKRwM3NDbGxsZg6dSqcnJx6XRwRERFZ3/T4AGRsPIYdhTVoadPDzcWxvrN7HG7++c9/4uzZs2hqaoKPT8cFgOrq6uDh4QG5XI6qqipER0dj8+bNPIGbiIjIBg0OlCNY6YYKdQt2nqzB9PgAsUuyqh4vKH799dcxduxYnDhxAjU1NaipqcHx48cxfvx4LF++HCUlJQgKCsKTTz7ZF/USERFRL0kkEkyP77h8S3b+WZGrsb4eh5sXX3wR//znPxETE2O8LzY2Fn//+9+Rnp6OsLAwvPnmm9i+fbtVCyUiIiLrmTa4I9xsOc5wg4qKCrS3t190f3t7u/EKxSEhIdBqtb2vjoiIiPrEpFg/OEslKKpuxKkax9oS3uNwc+WVV+LBBx9ETk6O8b6cnBwsWrQIV111FQAgLy8PUVGOt/qaiIjIUXi5uRgPz3S0qakeh5sPP/wQKpUKSUlJkMlkkMlkuOKKK6BSqfDhhx8CAORyOf7xj39YvVgiIiKyns6FxNn5VSJXYl09PhW807Fjx3D8+HEAQHx8POLj461aWF/hqeBEREQdjlZokLp8K9xcpMhdPMumt4T35Pv7si/iN2TIEAwZMuRyn05EREQiGxLkhSCFGyo1LdhVVGtcZGzvLivclJWV4ZtvvkFJSQlaW1tNHnvrrbesUhgRERH1LYlEgmmD/fH53lJsyT87cMPNpk2bcMMNNyA6OhrHjh3D8OHDUVxcDEEQMGbMmL6okYiIiPrI9PiOcJN9vAqLMUzscqyixwuK09PT8cwzzyAvLw9ubm748ssvUVpaimnTpuHWW2/tixqJiIioj0yK84OTVIKTZxtRWtskdjlW0eNwc/ToUaSlpQEAnJ2d0dzcDLlcjldffRVvvPGG1QskIiKivqNwc0FSROeWcMfYNdXjcOPp6WlcZxMcHIzCwkLjY9XV1darjIiIiPrF5Dg/AMCe4jqRK7GOHq+5mTBhArZt24ahQ4fi2muvxdNPP428vDysX78eEyZM6IsaiYiIqA+NDPcGAOSVq8UtxEp6HG7eeustNDQ0AACWLFmChoYGfP7554iLi+NOKSIiIjs0IlQJACiqboSmpQ0KNxeRK+qdHoeb6Oho4589PT3x73//26oFERERUf9SeboizMcdZXXNOFSuxsQYP7FL6pUer7mJjo5GTU3NRffX19ebBB8iIiKyH52jN3ll9j811eNwU1xcDL1ef9H9Op0O5eXlVimKiIiI+teIsI5wc9AB1t1YPC31zTffGP+clZUFpVJp/Fmv12PTpk0YNGiQVYsjIiKi/pEY6g0AODSQws3cuXMBdFyqed68eSaPubi4YNCgQTwJnIiIyE51TkudqmmCuqkNSg/7XVRscbgxGAwAgKioKOzZswd+fva92IiIiIjOU3q4INLXA6dqmpBXrjZe+8Ye9XjNTVFREYMNERGRAxoe2rnupl7cQnrJopGbd955x+IXfOyxxy67GCIiIhJPYqgS3x+ssPt1NxaFm3/+858WvZhEImG4ISIislPGHVN2vh3conBTVFTU13UQERGRyDqnpcrqmlHb2AqVp6vIFV2eHq+5uZAgCBAEwVq1EBERkYgUbi6I9vMEYN/nTF1WuPn4448xYsQIuLu7w93dHYmJifjkk0+sXRsRERH1s87RG3ted3NZB2e+9NJLeOSRRzBp0iQAwLZt27Bw4UJUV1fjySeftHqRRERE1D8Sw5T45sBpHCyrF7uUy9bjcPPuu+8iMzMTaWlpxvtuuOEGJCQk4JVXXmG4ISIismOOcMZUj6elKioqMHHixIvunzhxIioqKqxSFBEREYkjIVQJiQQ4rW5BdYNO7HIuS4/DTWxsLP73v/9ddP/nn3+OuLg4qxRFRERE4pDLnBHjLwdgv4uKezwttWTJEtx222347bffjGtutm/fjk2bNnUZeoiIiMi+jAhVoqCqAXllalwZHyB2OT1m8cjNoUOHAAA333wzdu3aBT8/P3z11Vf46quv4Ofnh927d+PGG2/ss0KJiIiof3Suu7HXi/lZPHKTmJiIsWPH4v7778ftt9+OTz/9tC/rIiIiIpEknrtScZ6dnjFl8cjNli1bkJCQgKeffhrBwcGYP38+tm7d2pe1ERERkQiGhSgglQBnNDqc1drfomKLw82UKVPw0UcfoaKiAu+++y6Kioowbdo0DB48GG+88QYqKyv7sk4iIiLqJx6uzgj1cQcAnDzbIHI1Pdfj3VKenp645557sGXLFhw/fhy33nor3nvvPUREROCGG27oixqJiIionw3y7TiGobimUeRKeq5XZ0vFxsbi+eefx4svvggvLy98//331qqLiIiIRBR17oypouomkSvpucsON7/99hvmz5+PoKAgPPvss7jpppuwffv2Hr1GZmYmEhMToVAooFAokJycjI0bN1r03LVr10IikWDu3LmXUT0RERGZ0xluiqvtb+SmR9e5OX36NFatWoVVq1ahoKAAEydOxDvvvIM//elP8PT07PGbh4WFYenSpYiLi4MgCFi9ejXmzJmDnJwcJCQkdPu84uJiPPPMM5gyZUqP35OIiIgubZCf/U5LWRxuUlNT8csvv8DPzw9paWm49957ER8f36s3nz17tsnPr732GjIzM7Fz585uw41er8ddd92FJUuWYOvWraivr+9VDURERHSxqAvW3BgMAqRSicgVWc7icOPi4oIvvvgC119/PZycnKxeiF6vx7p169DY2Ijk5ORu27366qsICAjAfffdZ9FWdJ1OB53u/DY2jUZjlXqJiIgcWZiPO5ylErS0GXBG24JgpbvYJVnM4nDzzTff9EkBeXl5SE5ORktLC+RyOTZs2IBhw4Z12Xbbtm348MMPkZuba/HrZ2RkYMmSJVaqloiIaGBwdpIiXOWBoupGFFU32lW46dVuKWuIj49Hbm4udu3ahUWLFmHevHk4cuTIRe20Wi3uvvtufPDBB/Dz87P49dPT06FWq4230tJSa5ZPRETksAb5egAAiu1sx1SPD860NldXV8TGxgIAkpKSsGfPHixfvhwrVqwwaVdYWIji4mKTdToGgwEA4OzsjPz8fMTExFz0+jKZDDKZrA97QERE5JgG+XkC+WftblGx6OHmjwwGg8kamU5DhgxBXl6eyX0vvvgitFotli9fjvDw8P4qkYiIaEA4f60bhhuLpaenIzU1FREREdBqtVizZg2ys7ORlZUFAEhLS0NoaCgyMjLg5uaG4cOHmzzf29sbAC66n4iIiHrPeJVihhvLVVVVIS0tDRUVFVAqlUhMTERWVhZmzpwJACgpKYFUKvqyICIiogGpc+TmVG2TXW0HlwiCIIhdRH/SaDRQKpVQq9VQKBRil0NERGSz9AYBQ1/6Ea16A7Y9dyXCfDxEq6Un398cFiEiIqIuOUklCFd1bAG3px1TDDdERETULeOiYjvaMcVwQ0RERN2yx0XFDDdERETUrUF2eDo4ww0RERF1i9NSRERE5FA6R25Ka5vQrjeIXI1lGG6IiIioW8EKN8icpWjTCzhd3yJ2ORZhuCEiIqJuSaUSRJ47QNNepqYYboiIiMgse9sxxXBDREREZtnbAZoMN0RERGSWcTs4p6WIiIjIEXBaioiIiBxK57RUaV0z2uxgOzjDDREREZkVqJDB3cUJeoOAsrpmscu5JIYbIiIiMksiuWA7eHWDyNVcGsMNERERXdL5HVNNIldyaQw3REREdEn2dIAmww0RERFdUqSqY1qqrI4jN0REROQAApVuAIAzGp3IlVwaww0RERFdUqBXZ7ix/cMzGW6IiIjokoLOjdzUNLZC164XuRrzGG6IiIjoknw8XODq1BEbzmpte2qK4YaIiIguSSKRIEAhA2D7U1MMN0RERGSRIEXH1FSlmiM3RERE5ADO75jiyA0RERE5AHvZMcVwQ0RERBYJUnasualkuCEiIiJHEKjgyA0RERE5kCCFfVylmOGGiIiILBJo3C3VAkEQRK6meww3REREZJHOqxQ3t+mhaWkXuZruMdwQERGRRdxcnKB0dwEAVNnwuhuGGyIiIrJYoML2d0wx3BAREZHFLlx3Y6sYboiIiMhinTumqmz48EyGGyIiIrIYR26IiIjIoXSeL8U1N0REROQQjNNSDDdERETkCLhbioiIiBxK58jNWa0O7XqDyNV0jeGGiIiILOYrl8FJKoFBAGoaW8Uup0sMN0RERGQxJ6kE/vJzU1M2umOK4YaIiIh6xNZ3TDHcEBERUY8EnVtUbKs7phhuiIiIqEc6FxVz5KYLmZmZSExMhEKhgEKhQHJyMjZu3Nht+/Xr1+OKK66At7c3PD09MWrUKHzyySf9WDEREREFGK9SbJtHMDiL+eZhYWFYunQp4uLiIAgCVq9ejTlz5iAnJwcJCQkXtVepVHjhhRcwZMgQuLq64rvvvsM999yDgIAAXH311SL0gIiIaODpHLk5Y6MjNxJBEASxi7iQSqXCsmXLcN9991nUfsyYMbjuuuvw17/+1aL2Go0GSqUSarUaCoWiN6USERENSNsLqnHXyl2IC5Dj56em9ct79uT722bW3Oj1eqxduxaNjY1ITk6+ZHtBELBp0ybk5+dj6tSp3bbT6XTQaDQmNyIiIrp8tn6VYlGnpQAgLy8PycnJaGlpgVwux4YNGzBs2LBu26vVaoSGhkKn08HJyQnvv/8+Zs6c2W37jIwMLFmypC9KJyIiGpA6TwbXtrSjqbUdHq6ixwkToo/cxMfHIzc3F7t27cKiRYswb948HDlypNv2Xl5eyM3NxZ49e/Daa6/hqaeeQnZ2drft09PToVarjbfS0tI+6AUREdHA4eXmAk9XJwDAGY3tLSq2uTU3KSkpiImJwYoVKyxqf//996O0tBRZWVkWteeaGyIiot676u/ZOFndiM8WTEByjG+fv59drrnpZDAYoNNZngJ72p6IiIh6L9CGd0yJOkmWnp6O1NRUREREQKvVYs2aNcjOzjaOwqSlpSE0NBQZGRkAOtbPXHHFFYiJiYFOp8MPP/yATz75BJmZmWJ2g4iIaMAJUjLcdKmqqgppaWmoqKiAUqlEYmIisrKyjAuES0pKIJWeH1xqbGzEQw89hLKyMri7u2PIkCH49NNPcdttt4nVBSIiogEpwIZ3TNncmpu+xjU3REREvfef7UVY8u0RXDsiCO/fldTn72fXa26IiIjI9p2/SrHtrXtluCEiIqIeO3++lO1NSzHcEBERUY91Liiu0rbAYLCtFS4MN0RERNRjAV4dC4rb9ALqmlpFrsYUww0RERH1mIuTFH5yVwC2t2OK4YaIiIgui5+8Y/SmtpEjN0REROQAfDw6Rm4YboiIiMgh+Hi6AADqm9pErsQUww0RERFdFm+O3BAREZEjUZ0LN/XcLUVERESOwNujY1qqjtNSRERE5AhUnh0jN7zODRERETmEzt1SDDdERETkEIzTUo2cliIiIiIHwGkpIiIiciidW8GbWvVoadOLXM15DDdERER0WRRuznCSSgDY1oX8GG6IiIjoskgkEvgYt4PbztQUww0RERFdts6pqTobukoxww0RERFdNpVxOzinpYiIiMgBdG4Hr+W0FBERETmCzu3g9ZyWIiIiIkfgzWkpIiIiciQqT+6WIiIiIgfibYPnSzHcEBER0WXz4VZwIiIiciTnp6W45oaIiIgcAC/iR0RERA6l8yJ+Wl072vQGkavpwHBDREREl03h7gJJx9mZNnN4JsMNERERXTYnqQRKd9vaDs5wQ0RERL2isrF1Nww3RERE1Cud50tx5IaIiIgcQuf5UrayHZzhhoiIiHrF1q5SzHBDREREveLTOS3FNTdERETkCHw4LUVERESOxNbOl2K4ISIiol7x4ZobIiIiciTGNTecliIiIiJHcH7NDUduiIiIyAF0Tkupm9ugNwgiV8NwQ0RERL3UeYViQegIOGJjuCEiIqJecXGSwsvNGYBtTE0x3BAREVGv2dJ2cFHDTWZmJhITE6FQKKBQKJCcnIyNGzd22/6DDz7AlClT4OPjAx8fH6SkpGD37t39WDERERF1xZYu5CdquAkLC8PSpUuxb98+7N27F1dddRXmzJmDw4cPd9k+Ozsbd9xxBzZv3owdO3YgPDwcs2bNQnl5eT9XTkRERBfysaGTwSWCIIi/rPkCKpUKy5Ytw3333XfJtnq9Hj4+PvjXv/6FtLQ0i15fo9FAqVRCrVZDoVD0tlwiIiIC8OTnudiQU4701CF4cFqM1V+/J9/fzlZ/98uk1+uxbt06NDY2Ijk52aLnNDU1oa2tDSqVqts2Op0OOp3O+LNGo+l1rURERGTq/FWKB/i0FADk5eVBLpdDJpNh4cKF2LBhA4YNG2bRc5977jmEhIQgJSWl2zYZGRlQKpXGW3h4uLVKJyIionNs6WRw0cNNfHw8cnNzsWvXLixatAjz5s3DkSNHLvm8pUuXYu3atdiwYQPc3Ny6bZeeng61Wm28lZaWWrN8IiIigm1dpVj0aSlXV1fExsYCAJKSkrBnzx4sX74cK1as6PY5f//737F06VL88ssvSExMNPv6MpkMMpnMqjUTERGRKVs6PFP0cPNHBoPBZI3MH7355pt47bXXkJWVhSuuuKIfKyMiIqLu2NLhmaKGm/T0dKSmpiIiIgJarRZr1qxBdnY2srKyAABpaWkIDQ1FRkYGAOCNN97A4sWLsWbNGgwaNAiVlZUAALlcDrlcLlo/iIiIBrrOaan6gT5yU1VVhbS0NFRUVECpVCIxMRFZWVmYOXMmAKCkpARS6fllQZmZmWhtbcUtt9xi8jovv/wyXnnllf4snYiIiC5w4W4pQRAgkUhEq0XUcPPhhx+afTw7O9vk5+Li4r4rhoiIiC5b5+GZeoMATUs7lO4uotUi+m4pIiIisn9uLk7wcHUCIP52cIYbIiIisgpb2THFcENERERW4ePZMRVVL/KOKYYbIiIisorOkZtaTksRERGRI+C0FBERETmU8xfyY7ghIiIiB3D+fCmuuSEiIiIHYJyW4pobIiIicgTenJYiIiIiR6LydIWTVAJBELcOmzsVnIiIiOzTpBg/FLyWKuq5UgDDDREREVmJVCpuqOnEaSkiIiJyKAw3RERE5FAYboiIiMihMNwQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNwQERGRQ2G4ISIiIocy4E4FFwQBAKDRaESuhIiIiCzV+b3d+T1uzoALN1qtFgAQHh4uciVERETUU1qtFkql0mwbiWBJBHIgBoMBp0+fhpeXFyQSiVVfW6PRIDw8HKWlpVAoFFZ9bVs2UPsNsO8Dse8Dtd8A+z4Q+25L/RYEAVqtFiEhIZBKza+qGXAjN1KpFGFhYX36HgqFQvRfAjEM1H4D7PtA7PtA7TfAvg/EvttKvy81YtOJC4qJiIjIoTDcEBERkUNhuLEimUyGl19+GTKZTOxS+tVA7TfAvg/Evg/UfgPs+0Dsu732e8AtKCYiIiLHxpEbIiIicigMN0RERORQGG6IiIjIoTDcEBERkUNhuLGS9957D4MGDYKbmxvGjx+P3bt3i11Sr7zyyiuQSCQmtyFDhhgfb2lpwcMPPwxfX1/I5XLcfPPNOHPmjMlrlJSU4LrrroOHhwcCAgLw7LPPor29vb+7ckm//fYbZs+ejZCQEEgkEnz11VcmjwuCgMWLFyM4OBju7u5ISUnBiRMnTNrU1tbirrvugkKhgLe3N+677z40NDSYtDl48CCmTJkCNzc3hIeH48033+zrrl3Spfo+f/78i34PrrnmGpM29tj3jIwMjB07Fl5eXggICMDcuXORn59v0sZav+PZ2dkYM2YMZDIZYmNjsWrVqr7uXrcs6ff06dMv+swXLlxo0sbe+g0AmZmZSExMNF6MLjk5GRs3bjQ+7oifd6dL9d0hP3OBem3t2rWCq6ur8NFHHwmHDx8WFixYIHh7ewtnzpwRu7TL9vLLLwsJCQlCRUWF8Xb27Fnj4wsXLhTCw8OFTZs2CXv37hUmTJggTJw40fh4e3u7MHz4cCElJUXIyckRfvjhB8HPz09IT08Xoztm/fDDD8ILL7wgrF+/XgAgbNiwweTxpUuXCkqlUvjqq6+EAwcOCDfccIMQFRUlNDc3G9tcc801wsiRI4WdO3cKW7duFWJjY4U77rjD+LharRYCAwOFu+66Szh06JDw2WefCe7u7sKKFSv6q5tdulTf582bJ1xzzTUmvwe1tbUmbeyx71dffbXwn//8Rzh06JCQm5srXHvttUJERITQ0NBgbGON3/GTJ08KHh4ewlNPPSUcOXJEePfddwUnJyfhxx9/7Nf+drKk39OmTRMWLFhg8pmr1Wrj4/bYb0EQhG+++Ub4/vvvhePHjwv5+fnC888/L7i4uAiHDh0SBMExP+9Ol+q7I37mDDdWMG7cOOHhhx82/qzX64WQkBAhIyNDxKp65+WXXxZGjhzZ5WP19fWCi4uLsG7dOuN9R48eFQAIO3bsEASh40tTKpUKlZWVxjaZmZmCQqEQdDpdn9beG3/8gjcYDEJQUJCwbNky43319fWCTCYTPvvsM0EQBOHIkSMCAGHPnj3GNhs3bhQkEolQXl4uCIIgvP/++4KPj49J35977jkhPj6+j3tkue7CzZw5c7p9jqP0vaqqSgAgbNmyRRAE6/2O/9///Z+QkJBg8l633XabcPXVV/d1lyzyx34LQscX3eOPP97tcxyh3518fHyElStXDpjP+0KdfRcEx/zMOS3VS62trdi3bx9SUlKM90mlUqSkpGDHjh0iVtZ7J06cQEhICKKjo3HXXXehpKQEALBv3z60tbWZ9HnIkCGIiIgw9nnHjh0YMWIEAgMDjW2uvvpqaDQaHD58uH870gtFRUWorKw06atSqcT48eNN+urt7Y0rrrjC2CYlJQVSqRS7du0ytpk6dSpcXV2Nba6++mrk5+ejrq6un3pzebKzsxEQEID4+HgsWrQINTU1xsccpe9qtRoAoFKpAFjvd3zHjh0mr9HZxlb+bfhjvzv997//hZ+fH4YPH4709HQ0NTUZH3OEfuv1eqxduxaNjY1ITk4eMJ83cHHfOznaZz7gDs60turqauj1epMPHQACAwNx7NgxkarqvfHjx2PVqlWIj49HRUUFlixZgilTpuDQoUOorKyEq6srvL29TZ4TGBiIyspKAEBlZWWXfyedj9mLzlq76suFfQ0ICDB53NnZGSqVyqRNVFTURa/R+ZiPj0+f1N9b11xzDW666SZERUWhsLAQzz//PFJTU7Fjxw44OTk5RN8NBgOeeOIJTJo0CcOHDzfWZY3f8e7aaDQaNDc3w93dvS+6ZJGu+g0Ad955JyIjIxESEoKDBw/iueeeQ35+PtavXw/Avvudl5eH5ORktLS0QC6XY8OGDRg2bBhyc3Md/vPuru+AY37mDDfUpdTUVOOfExMTMX78eERGRuJ///ufqP8Dpf51++23G/88YsQIJCYmIiYmBtnZ2ZgxY4aIlVnPww8/jEOHDmHbtm1il9Kvuuv3Aw88YPzziBEjEBwcjBkzZqCwsBAxMTH9XaZVxcfHIzc3F2q1Gl988QXmzZuHLVu2iF1Wv+iu78OGDXPIz5zTUr3k5+cHJyeni1bVnzlzBkFBQSJVZX3e3t4YPHgwCgoKEBQUhNbWVtTX15u0ubDPQUFBXf6ddD5mLzprNff5BgUFoaqqyuTx9vZ21NbWOtzfR3R0NPz8/FBQUADA/vv+yCOP4LvvvsPmzZsRFhZmvN9av+PdtVEoFKL+n4Tu+t2V8ePHA4DJZ26v/XZ1dUVsbCySkpKQkZGBkSNHYvny5Q7/eQPd970rjvCZM9z0kqurK5KSkrBp0ybjfQaDAZs2bTKZz7R3DQ0NKCwsRHBwMJKSkuDi4mLS5/z8fJSUlBj7nJycjLy8PJMvvp9//hkKhcI4FGoPoqKiEBQUZNJXjUaDXbt2mfS1vr4e+/btM7b59ddfYTAYjP9IJCcn47fffkNbW5uxzc8//4z4+HjRp2V6oqysDDU1NQgODgZgv30XBAGPPPIINmzYgF9//fWiaTNr/Y4nJyebvEZnG7H+bbhUv7uSm5sLACafub31uzsGgwE6nc5hP29zOvveFYf4zEVZxuxg1q5dK8hkMmHVqlXCkSNHhAceeEDw9vY2WVlub55++mkhOztbKCoqErZv3y6kpKQIfn5+QlVVlSAIHdsmIyIihF9//VXYu3evkJycLCQnJxuf37l1cNasWUJubq7w448/Cv7+/ja5FVyr1Qo5OTlCTk6OAEB46623hJycHOHUqVOCIHRsBff29ha+/vpr4eDBg8KcOXO63Ao+evRoYdeuXcK2bduEuLg4k+3Q9fX1QmBgoHD33XcLhw4dEtauXSt4eHiIvhXcXN+1Wq3wzDPPCDt27BCKioqEX375RRgzZowQFxcntLS0GF/DHvu+aNEiQalUCtnZ2SbbX5uamoxtrPE73rk99tlnnxWOHj0qvPfee6Juj71UvwsKCoRXX31V2Lt3r1BUVCR8/fXXQnR0tDB16lTja9hjvwVBEP7yl78IW7ZsEYqKioSDBw8Kf/nLXwSJRCL89NNPgiA45ufdyVzfHfUzZ7ixknfffVeIiIgQXF1dhXHjxgk7d+4Uu6Reue2224Tg4GDB1dVVCA0NFW677TahoKDA+Hhzc7Pw0EMPCT4+PoKHh4dw4403ChUVFSavUVxcLKSmpgru7u6Cn5+f8PTTTwttbW393ZVL2rx5swDgotu8efMEQejYDv7SSy8JgYGBgkwmE2bMmCHk5+ebvEZNTY1wxx13CHK5XFAoFMI999wjaLVakzYHDhwQJk+eLMhkMiE0NFRYunRpf3WxW+b63tTUJMyaNUvw9/cXXFxchMjISGHBggUXhXZ77HtXfQYg/Oc//zG2sdbv+ObNm4VRo0YJrq6uQnR0tMl79LdL9bukpESYOnWqoFKpBJlMJsTGxgrPPvusyTVPBMH++i0IgnDvvfcKkZGRgqurq+Dv7y/MmDHDGGwEwTE/707m+u6on7lEEASh/8aJiIiIiPoW19wQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FAYboiIiMihMNwQkejmz5+PuXPn9vv7rlq1ChKJBBKJBE888YRFz5k/f77xOV999VWf1kdEl8dZ7AKIyLFJJBKzj7/88stYvnw5xLpYukKhQH5+Pjw9PS1qv3z5cixdutR4qCAR2R6GGyLqUxUVFcY/f/7551i8eDHy8/ON98nlcsjlcjFKA9ARvoKCgixur1QqoVQq+7AiIuotTksRUZ8KCgoy3pRKpTFMdN7kcvlF01LTp0/Ho48+iieeeAI+Pj4IDAzEBx98gMbGRtxzzz3w8vJCbGwsNm7caPJehw4dQmpqKuRyOQIDA3H33Xejurq6xzW///77iIuLg5ubGwIDA3HLLbf09q+BiPoRww0R2aTVq1fDz88Pu3fvxqOPPopFixbh1ltvxcSJE7F//37MmjULd999N5qamgAA9fX1uOqqqzB69Gjs3bsXP/74I86cOYM//elPPXrfvXv34rHHHsOrr76K/Px8/Pjjj5g6dWpfdJGI+ginpYjIJo0cORIvvvgiACA9PR1Lly6Fn58fFixYAABYvHgxMjMzcfDgQUyYMAH/+te/MHr0aLz++uvG1/joo48QHh6O48ePY/DgwRa9b0lJCTw9PXH99dfDy8sLkZGRGD16tPU7SER9hiM3RGSTEhMTjX92cnKCr68vRowYYbwvMDAQAFBVVQUAOHDgADZv3mxcwyOXyzFkyBAAQGFhocXvO3PmTERGRiI6Ohp33303/vvf/xpHh4jIPjDcEJFNcnFxMflZIpGY3Ne5C8tgMAAAGhoaMHv2bOTm5prcTpw40aNpJS8vL+zfvx+fffYZgoODsXjxYowcORL19fW97xQR9QtOSxGRQxgzZgy+/PJLDBo0CM7OvfunzdnZGSkpKUhJScHLL78Mb29v/Prrr7jpppusVC0R9SWO3BCRQ3j44YdRW1uLO+64A3v27EFhYSGysrJwzz33QK/XW/w63333Hd555x3k5ubi1KlT+Pjjj2EwGBAfH9+H1RORNTHcEJFDCAkJwfbt26HX6zFr1iyMGDECTzzxBLy9vSGVWv5Pnbe3N9avX4+rrroKQ4cOxb///W989tlnSEhI6MPqiciaJIJYlwUlIhLZqlWr8MQTT1zWehqJRIINGzaIcmwEEZnHkRsiGtDUajXkcjmee+45i9ovXLhQ1CsqE9GlceSGiAYsrVaLM2fOAOiYjvLz87vkc6qqqqDRaAAAwcHBFp9JRUT9h+GGiIiIHAqnpYiIiMihMNwQERGRQ2G4ISIiIofCcENEREQOheGGiIiIHArDDRERETkUhhsiIiJyKAw3RERE5FD+P1bdBYJlA0RMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create simulation\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "sim = pybamm.Simulation(model, solver=solver)\n",
    "\n",
    "# solve\n",
    "solution = sim.solve([0, 3600])\n",
    "\n",
    "# Post-processing\n",
    "t_evals = np.linspace(0, 3600, 100)\n",
    "voltage = solution[\"Terminal voltage [V]\"](t_evals)\n",
    "\n",
    "# Plot the results\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(t_evals, voltage)\n",
    "plt.xlabel(\"Time [s]\")\n",
    "plt.ylabel(\"Voltage [V]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we can instrument this script to time the different parts. In this case, the entire pipeline is performed during the `solve` call, but helpfully the `pybamm.Simulation` class will cache the discretisation and solver setup so that they only need to be performed once. So subsequent calls to `solve` will be faster than the first call and will only include the solver solve step (step 4). So we can time the first call to `solve` and then time subsequent calls to `solve` to see the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First solve took 0.411 seconds\n",
      "Second solve took 0.026 seconds\n",
      "Third solve took 0.026 seconds\n",
      "Post-processing took 0.002 seconds\n"
     ]
    }
   ],
   "source": [
    "model = pybamm.lithium_ion.DFN()\n",
    "solver = pybamm.IDAKLUSolver()\n",
    "sim = pybamm.Simulation(model, solver=solver)\n",
    "\n",
    "# solve\n",
    "start = time.perf_counter()\n",
    "solution = sim.solve([0, 3600])\n",
    "end = time.perf_counter()\n",
    "print(f\"First solve took {end - start:.3f} seconds\")\n",
    "\n",
    "# solve again\n",
    "start = time.perf_counter()\n",
    "solution = sim.solve([0, 3600])\n",
    "end = time.perf_counter()\n",
    "print(f\"Second solve took {end - start:.3f} seconds\")\n",
    "\n",
    "# solve again\n",
    "start = time.perf_counter()\n",
    "solution = sim.solve([0, 3600])\n",
    "end = time.perf_counter()\n",
    "print(f\"Third solve took {end - start:.3f} seconds\")\n",
    "\n",
    "# Post-processing\n",
    "start = time.perf_counter()\n",
    "t_evals = np.linspace(0, 3600, 100)\n",
    "voltage = solution[\"Terminal voltage [V]\"](t_evals)\n",
    "end = time.perf_counter()\n",
    "print(f\"Post-processing took {end - start:.3f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this output you can see that the first call to `solve` takes roughly the same amount of time as the entire pipeline, but subsequent calls to `solve` are much faster because the discretisation and solver setup have already been performed. \n",
    "\n",
    "## Conclusion\n",
    "\n",
    "The main takeaway from this is that the simulation object is a heavy object that caches many of the intermediate steps of the pipeline, so if you are running many simulations with the same model and discretisation, you should use the same simulation object to avoid the overhead of setting up the simulation each time. For the particular case of varying the parameters of the model, you can use input parameters to keep certain model parameters as symbols rather than concrete values during the 1st step of the pipeline, thus allowing you to re-use a particular discretisation and solver setup. This will be discussed in more depth in the next notebook."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
